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ABSTRACT 


Signal processing methods for signals sampled at different rates are inves¬ 
tigated and applied to the problem of signal and image reconstruction or super¬ 
resolution reconstruction. The problem is approached from the viewpoint of linear 
mean-square estimation theory and multirate signal processing for one- and two- 
dimensional signals. A new look is taken at multirate system theory in one and two 
dimensions which provides the framework for these methodologies. A careful analysis 
of linear optimal hltering for problems involving different input and output sampling 
rates is conducted. This results in the development of index mapping techniques that 
simplify the formulation of Wiener-Hopf equations whose solution determine the op¬ 
timal hlters. The required hlters exhibit periodicity in both one and two dimensions, 
due to the difference in sampling rates. The reconstruction algorithms developed are 
applied to one- and two-dimensional reconstruction problems. 
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EXECUTIVE SUMMARY 


As physical and manufacturing limitations are reached in state-of-the-art im¬ 
age acquisition systems, there is increased motivation to improve the resolution of 
imagery through signal processing methods. High-resolution (HR) imagery is desir¬ 
able because it can offer more detail about the object associated with the imagery. 
The “extra” information is of critical importance in many applications. For exam¬ 
ple, HR reconnaissance images can provide intelligence analysts, greater information 
about a military target, including its capabilities, operability and vulnerabilities, and 
increase analysts’ conhdence in such assessments. Likewise, HR medical images can 
be crucial to a physician in making a proper diagnosis or developing a suitable treat¬ 
ment regimen. 

Super-resolution (SR) image reconstruction is an approach to this problem, 
and this area of research encompasses those signal processing techniques that use 
multiple low-resolution (LR) images to form a HR image of some related object. In 
this work, a super-resolution image reconstruction approach is proposed from the 
viewpoint of estimation and multirate signal processing for two-dimensional signals 
or images. 

Multirate signal processing theory deals with the analysis of a system com¬ 
prised of multiple signals at different sampling rates and is fundamental to this re¬ 
search. An example of such a system is a sensor network that collects and processes 
data from various sensors, where the information from each sensor might be collected 
at a different rate. In developing this theory, a number of relationships between sig¬ 
nals in a multirate system are identihed. The critical hnding is that all of the signals 
in a multirate system can be referred to a single “universal” rate for that system; 
therefore, many of the results of standard signal processing theory can be adapted to 
multirate systems through this observation. 


xvii 



The multirate theory developed here is applied to signal estimation, where one 
signal is estimated from some other related signal or signals. The desired signal may 
be corrupted by distortion or interference and is usually unobservable (at least at 
the moment when the estimate is desired). A typical signal estimation application is 
the recovery of a transmitted signal from a received signal that has been subject to 
distortion and is corrupted by noise. 

SR image reconstruction can be viewed as a problem in signal estimation, 
where a related LR signal or signals is used to estimate an underlying HR signal. 
From this perspective, the observation signal or signals, and desired signal form a 
multirate system. This motivates the application of the theory of multirate systems 
to signal estimation and the resultant extension of single-rate signal estimation theory 
to the multirate case. 

The particular branch of estimation theory applied in this work is optimal 
hltering, where the error in estimation is minimized by using a weighted set of the 
LR observation images to hlter and estimate the HR image. The weights used in this 
linear estimate are called hlter coefficients and application of this theory results in 
a set of equations that are solved to obtain these coefficients known as the Wiener- 
Hopf (WH) equations. In this research, the multirate WH equations are developed 
and shown to have a periodically time-dependent solution. Additionally, the concept 
of index mapping, an extension of the multirate theory, is developed to determine the 
required regions of the LR images required for estimation. 

A new methodology is developed and presented, by application and extension 
of the results of multirate and optimal estimation theory to the problem of SR image 
reconstruction. This new method is applied to a set of LR images, and the resultant 
HR image is compared with results from standard interpolation methods. In every 
case, this method performed better than the standard methods. 
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I. INTRODUCTION 


As physical and manufacturing limitations are reached in state-of-the-art im¬ 
age acquisition systems, there is increased motivation to improve the resolution of 
imagery through signal processing methods. Improvements in this area have signih- 
cant commercial and military application, and in this work a super-resolution image 
reconstruction approach is proposed from the viewpoint of estimation and multirate 
signal processing for two-dimensional signals. 

A. PROBLEM STATEMENT/MOTIVATION 

Super-resolution (SR) imaging has recently become an area of great interest 
in the image processing research community (see Section I.B.2). The ability to form 
a high-resolution (HR) image from a collection of subsampled images has a broad 
range of applications and has largely been motivated by physical and production 
limitations on existing image acquisition systems and the marginal costs associated 
with increased spatial resolution. Figure 1.1 depicts the SR concept where a collection 
of low-resolution (LR) images of a scene are superimposed on a HR grid, available 
for subsequent HR image reconstruction. 

In this work, we propose a stochastic multirate approach to this problem, 
adapting and extending the work in [Ref. 6, 7, 8, 9] to one- and two-dimensional 
signals. The earlier work has focused on information fusion applications, i.e., on the 
combination of observations from multiple sensors to perform tracking, surveillance, 
classihcation or some other task. This work extends these concepts to reconstruction 
of one-dimensional signals and SR image reconstruction. 
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Figure 1.1. Super-resolution imaging concept, (After [Ref. 1]). 


B. PREVIOUS WORK 

1. Stochastic Multirate Signal Processing 

Research in the area of stochastic multirate signal processing has been lim¬ 
ited to a handful of investigators whose work has focused mainly on second moment 
analysis of stochastic systems, from both temporal and spectral points of view, and 
optimal estimation theory, including both Kalman and Weiner filtering theory. 

Vaidyanathan et al. [Ref. 10, 11, 12] investigate how the statistical properties 
of stochastic signals are altered through multirate systems. In [Ref. 10], several facts 
and theorems are presented regarding the statistical behavior of signals as they are 
passed through decimators, interpolators, modulators, and more complicated inter- 
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connections. For example, the necessary and sufficient condition for the output of 
an L-fold interpolation hlter to be wide-sense stationary (WSS), given a WSS input, 
is that the L-fold decimation of the hlter coefficients results in no aliasing, i. e., the 
hlter must have an alias-free (L) support. Additionally, the authors illustrate an 
application of this theoretical analysis to a multirate adaptive hltering scheme for 
identihcation of band-limited channels. In [Ref. 11], this work is continued but ad¬ 
dressed using bifrequency maps and bispectra. These two-dimensional (2-D) Fourier 
transforms characterize all linear time-varying (LTV) systems and nonstationary ran¬ 
dom processes, respectively. In fact, by using these concepts, the previous results are 
simplihed and even generalized to handle the case of vector systems. Finally, in 
[Ref. 12], further analysis is conducted using bifrequency maps and bispectra, and a 
bifrequency characterization of lossless LTV systems is derived. 

Jahromi et al. [Ref. 13, 14, 15] consider methods to optimally estimate samples 
of a random signal based on observations made by multiple observers at diherent 
sampling rates (lower than the original rate). In particular, in [Ref. 13], the problem 
of fusing two low-rate sensors in the reconstruction of one high-resolution signal is 
considered when time delay of arrival (TDOA) is present. Using the “generalized 
cross-correlation” technique, the delay is estimated and then signal reconstruction is 
accomplished using perfect reconstruction synthesis filter bank theory. In [Ref. 14] 
and [Ref. 15], optimal least mean-square estimation is used to develop an estimate 
for samples of a high-rate signal. The estimator is a function of the power spectral 
density of the original random signal, which is obtained using a method for inductive 
inference of probability distribution referred to as the “maximum entropy principle” 
[Ref. 16]. 

Chen et al. [Ref. 17, 18, 19, 20] investigate use of the Kalman hlter and 
Weiner hlter in the reconstruction of a stochastic signal when only a noisy, downsam¬ 
pled version of the signal can be measured. In [Ref. 17], the use of the Kalman hlter 
is investigated for interpolating and estimating values of an autoregressive or moving 
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average stochastic signal when only a noisy, downsampled version of the signal can 
be measured. The signal reconstruction problem is converted into a state estima¬ 
tion problem for which the Kalman hlter is optimal. Some extensions are discussed, 
including the application of the Kalman reconstruction hlter in recovering missing 
speech packets in a packet switching network with packet interleaving. Simulation 
results are presented, which indicate that the multirate Kalman reconstruction hlters 
possess better reconstruction performance than a Wiener reconstruction hlter under 
comparable numerical complexity. In [Ref. 18], a multirate deconvolution hlter is pro¬ 
posed for signal reconstruction in multirate systems with channel noise. Both hlter 
bank and transmultiplexer architectures are used to demonstrate the design proce¬ 
dure. In [Ref. 19], a block state-space model is introduced where transmultiplexer 
systems unify the multirate signals and channel noise. In [Ref. 20], the optimal signal 
reconstruction problem is considered in transmultiplexer systems under channel noise 
from the viewpoint of Wiener-Hopf theory. A calculus of variation method and a 
spectral factorization technique are used to develop an appropriate separation hlter 
bank design. 

Scharf et al. [Ref. 21] introduce a least squares design methodology for hl- 
tering periodically correlated (PC) scalar time series. Since any PC time series can 
be represented as a WSS vector time series where each constituent subsequence is 
a decimated version of the original shifted in time, and vice versa, multirate hlter 
banks and equivalent polyphase realizations provide a natural representation for this 
bidirectional relationship. This relationship ahords means to develop a spectral rep¬ 
resentation for the PC time series and hence develop causal synthesis and causal 
whitening hlters for the PC scalar time series. These techniques are used to solve 
generalized linear minimum mean-square error (MMSE) hlter design problems for 
PC scalar time series. Note that this viewpoint can be extended to multirate systems 
where the correlation between observation sequences is periodically correlated. 
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Therrien et al. [Ref. 6, 22, 7, 8, 9, 23] develop theory and methodology 
required for employing optimal linear hltering in estimating an underlying signal 
from observation sequences at different sampling rates. The focus of these efforts is 
on information fusion, i.e., on the combination of observations from multiple sensors 
to perform tracking, surveillance, classihcation or some other task. In particular, 
[Ref. 6], [Ref. 22] and [Ref. 7] consider a simplihed problem where an underlying 
signal is estimated from two sequences, one observed at full rate and the other at 
half the rate. In [Ref. 8], least squares formulations are examined where the second 
sequence has an arbitrary sampling rate. In [Ref. 9], a general approach is suggested 
for any number of observation signals at arbitrary sampling rates. Finally, in [Ref. 
23], previous theory and methods are developed to consider the problem of HR signal 
and image reconstruction. This work forms the basis for the proposed research and 
represents an advance in the area of super-resolution image reconstruction. 

2. Super-Resolution Reconstruction/Imaging 

Generally, super-resolution (SR) image reconstruction refers to signal process¬ 
ing methods in which a high-resolution (HR) image is obtained from a set or ensemble 
of observed low-resolution (LR) images [Ref. 1]. If each observed LR image is sub¬ 
sampled (and aliased) and is translated by a different subpixel amount, this set of 
unique observation images can be used for reconstruction. Figure 1.1 demonstrates 
this conceptually. Both [Ref. 1] and [Ref. 2] provide general surveys of research to 
date regarding this topic, and the following major areas of research are identified: 
nonuniform interpolation, frequency domain, regularized SR reconstruction, projec¬ 
tion onto convex sets (POCS), maximum likelihood (ML) projection onto convex sets 
(ML-POCS) hybrid reconstruction, and other approaches [Ref. 1]. 

The most prevalent approaches in the literature are those based on nonuni¬ 
form interpolation. These approaches typically use a three-stage sequential process, 
comprised of registration, interpolation, and restoration. The registration step is a 
mapping of pixels from each LR image to a reference grid, which results in a HR grid 
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comprised of a set of nonuniformly spaced pixels. The interpolation step conforms 
these nonuniformly spaced pixels to a uniform sampling grid, which results in the 
upsampled HR image. Finally, restoration removes the effects of sensor distortion 
and noise. This scheme is depicted in Figure 1.2. Representative works include [Ref. 
24, 25, 26, 27]. 


LR Images 



HR Image 

y 


Figure 1.2. Typical model for nonuniform interpolation approach to SR, (From [Ref. 

2 ]). 


The frequency-domain approaches exploit the relationship between the discrete 
Fourier transforms (DFT) of the LR images and the continuous Fourier transform 
(CFT) of the desired HR image by using the information generated through relative 
motion between the LR images, the aliasing generated by downsampling relative to 
the desired HR image, and the assumption that the original HR image is bandlim- 
ited. A set of linear system equations are developed, and the continuous Fourier 
coefficients are found. The desired HR image is estimated from the CFT synthesis 
equation. Tsai and Huang [Ref. 28] were the hrst to introduce this method and were 
also the hrst researchers to address the problem of reconstructing a HR image from a 
set of translated LR images. Kim et al. [Ref. 29] extended this approach to include 
the presence of noise in the LR images using a recursive procedure based on weighted 
least squares theory. Kim and Su [Ref. 30] further extended this approach by consid- 
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ering noise and different blur distortions in the LR images. Vandewalle et al. [Ref. 
31] consider offset estimation using a subspace minimization method followed by a 
frequency-based reconstruction method based on the continuous and discrete Fourier 
series. 

The regularized SR reconstruction methods use regularization methods to solve 
the often ill-posed inverse problem introduced in the frequency-domain approaches. 
Typically, the ill-posed problems are a result of an insufficient number of LR images 
or ill-conditioned blur operators [Ref. 1]. Generally, two approaches have been con¬ 
sidered: deterministic and stochastic regularization. Deterministic approaches [Ref. 
32, 33, 34, 35] typically use constrained least squares methods (CLS) while stochastic 
approaches [Ref. 36, 37, 38] typically use maximum a posteriori (MAP) or maximum 
likelihood (ML) methods. 

ROCS methods are based on set theoretic estimation theory [Ref. 39]. Rather 
than using conventional estimation theory, the POCS formulations incorporate a pri¬ 
ori knowledge into the solution and yield a solution consistent with user-furnished 
constraints. Application of this method as applied to SR was introduced by Stark 
and Oskoui [Ref. 40] and extended by Tekalp et al. in [Ref. 41, 42], which takes 
into account the presence of both sensor blurring and observation noise, and suggests 
POCS as a new method for restoration of spatially-variant blurred images. 

ML-POCS hybrid reconstruction approaches estimate desired HR images by 
minimizing the ML or MAP cost functional while constraining the solution within 
certain closed convex sets in accordance with POCS methodology [Ref. 37]. 

There are a number of other areas that are considered in the literature, and 
some examples are presented here. One approach attempts to reconstruct a HR image 
from a single LR image and is referred to as improved dehnition image interpolation 
[Ref. 43]. Another area of study, referred to as iterative back-projection [Ref. 44, 45, 
46], uses tomographic projection methods to estimate a HR image. Researchers are 
also considering the SR problem when no relative subpixel motion exists between LR 
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images. By considering differently blurred LR images, motionless SR reconstruction 
can be demonstrated [Ref. 47, 48]. Milanfar et al. analyze the joint problem of 
image registration and HR reconstruction in the context of fundamental statistical 
performance limits. By using the Cramer-Rao bound, they demonstrate ability to 
bound estimator performance in terms of MSE, examining performance limits as 
they relate to such imaging system parameters as the downsampling factor, signal-to- 
noise ratio, and point spread function. Finally, researchers are considering adaptive 
hltering approaches to the SR problem, considering modihed recursive least squares 
(RLS), linear mean-square (LMS) and steepest descent methods [Ref. 49]. 

C. THESIS ORGANIZATION 

This manuscript is organized as follows. The current chapter is introductory 
and presents the motivation for this work, dehning the problem and outlining the 
approach used to solve it. Additionally, a review of the relevant literature is included, 
both in the area of stochastic multirate signal processing and super-resolution image 
reconstruction. 

The second chapter introduces various fundamental signal processing and 
mathematical concepts required for theoretic and application-related developments 
in future chapters. These include various signal taxonomies and representations, a 
review of relevant topics in second-moment analysis, and required number theory and 
linear algebra concepts. Further, this chapter, establishes notation and conventions 
for purposes of consistency throughout this work. 

In the third chapter, the theory of multirate systems is established. In this 
analysis, the relationships between a multirate system and its constituent signals are 
characterized, the system theory for multirate systems is developed, and the 



representation of discrete linear systems is presented from a system theoretic point 
of view. Finally, a linear algebraic approach is introduced to model various multirate 
operations for use in reconstruction applications. 

Chapter IV develops the concept of multirate signal estimation and is founda¬ 
tional in developing stochastic approaches to solving the signal reconstruction prob¬ 
lem. The optimal hltering problem is introduced in terms of the ordinary Wiener- 
Hopf equation and is then expanded, hrst to the single-channel, multirate estimation 
problem and then to the multi-channel, multirate problem. Also in this chapter, the 
relationship between samples in one signal domain to those in a different signal do¬ 
main (signals at different rates) is established through the concept of index mapping, 
which allows for a very general representation of the multirate Wiener-Hopf equations. 

Chapter V considers the problem of signal reconstruction in the one- and two- 
dimensions. In this chapter, the problem is stated for both cases, observation models 
are established, reconstruction approaches and algorithms are developed, and then 
the results of each algorithm are presented. 

Finally, Chapter VI provides conclusory remarks on the hndings of this re¬ 
search and establishes direction for future work related to this research. 
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II. PRELIMINARIES, CONVENTIONS, AND 

NOTATION 

In the development of approaches to signal and image reconstruction, a num¬ 
ber of fundamental concepts from the areas of signal processing and mathematics are 
required. In this chapter, a foundation is set in these areas upon which the theory of 
multirate signals and multirate estimation will be built. In doing so, we present the 
underlying concepts, but also emphasize required definitions, notations and conven¬ 
tions, in order to ensure consistency and accuracy, and to facilitate understanding. 

A. SIGNALS 

1. Etymology 

Etymologically speaking, the word signal is derived from the Latin signum, 
which can be rendered as “a sign, mark, or token;” or in a military sense, “a standard, 
banner, or ensign;” or “a physical representation of a person or thing, like a hgure, 
image, or statue [Ref. 50].” Generally, the Latin seems to imply that a signum is 
something that conveys information about or from someone or something else. The 
relevant modern dictionary dehnition of signal carries this idea further: “a detectable 
physical quantity or impulse by which messages or information can be transmitted 
[Ref. 51].” 

In the area of electrical engineering known as digital signal processing, a related 
but more helpful dehnition of a signal is a collection of information, usually a pattern 
of variation [Ref. 52], that describes some physical phenomenon. In other words, a 
signal conveys relevant information about some physical phenomena {signum). The 
variation in electrical voltage measured at the input of an electronic circuit, the 
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variation in acoustic pressure sensed by a microphone recording a musical concert, 
or the variation in light intensity captured by a camera recording a scene are all 
examples of signals treated in modern signal processing. 

2. Signal Definitions 

Throughout this presentation, various types of signals and sequences are in¬ 
troduced and analyzed. In this section, for the sake of clarity, the dehnition of such 
signals and sequences are established, as are the associated conventions and nota¬ 
tions. Let us begin with one-dimensional signals that are scalar-valued. We dehne 
these more precisely below. 

a. Deterministic Signals and Sequences 
A deterministic analog signal or simply an analog signal is dehned as 

follows. 

Definition 1. A deterministic analog signal, denoted by {a;(t)}, or when it is clear 
from context x{t), is a set of ordered measurements such that for every f G M, there 
exists a corresponding measurement m = x{t). If all such measurements are members 
of the extended real numbers^, then x(t) is said to be a real-valued (or real) analog 
signal. If the measurements are members of the complex numbers, then the signal is 
said to be a complex-valued (or complex) analog signal. 

An analog signal is frequently represented by a mathematical function, 
which may or may not be continuous. For example, the signal known as the unit-step, 
dehned by 

{ 1 t>0 

( 2 . 1 ) 

0 ( < 0 

is well known in signal processing, but the function representing it is not continuous 
(at t = 0). 

^The extended real numbers are defined as K = ]RU {—oo,oo}. 
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Although many signals are represented by functions defined on the real 
number line, our definition of a signal is not necessarily the same as the mathematical 
definition of a function. The set of analog signals commonly includes the unit impulse, 
which (strictly speaking) is not a function at all but a distribution or “generalized 
function,” described by a careful limiting process [Ref. 53, 54] to insure that the 
resulting entity satisfies certain conditions when it appears in an integral. 

Signals may have many other properties that provide for further char¬ 
acterization. One property of concern in this work is that of periodicity. A signal is 
said to be periodic if there exists a positive real number P such that 

x{t) = x{t + P) for all t. (2.2) 

The smallest such P is called the period. 

A deterministic sequence (or simply a sequence) is defined as follows. 

Definition 2. A deterministic sequence, denoted by {a:[n]}, or when clear from con¬ 
text x[n], is a countable set of ordered measurements such that for every n G Z, there 
exists a corresponding measurement m = x[n]. If all such measurements are members 
of the extended real numbers, then x[n\ is said to be a real-valued (or real) sequence. 
If the measurements are members of the complex numbers, then the sequence is said 
to be a complex-valued (or complex) sequence. 

A sequence x[n] is said to be periodic if there exists a positive integer 

N such that 

x[n] = x[n iV] for all n, (2.3) 

and the smallest such N is called the period. Note that not all sequences derived 
by sampling a periodic analog signal are periodic. For example, the analog signal 
x{t) = cos(27r/ot -|- 0) is periodic for any real number /o, while the sequence x[n] 
defined by x[n] = x{nTs) = cos{2nfouTg -|- 0) is periodic only if /qT* is a rational 
number. 
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Observe that both a signal and a sequence are dehned by an ordered 
set of measurements, but over a different domain (M or Z). Further, parentheses are 
used in the notation for an analog signal x{-) while square brackets are used for a 
sequence a;[-] (to indicate the discrete nature of its domain). The variable t or n is 
frequently used to represent time, although the units of “time” need to be specified 
in any real-world problem. In the case of a sequence, n is just an index variable used 
to order the measurements, and there is need in signal processing to dehne what will 
be called a deterministic discrete-domain signal or simply discrete-domain signal . 

Definition 3. A deterministic discrete-domain signal, denoted by {xT(t)}, or when 
it is clear from context XT{t), is a set of ordered measurements such that for every 
t G Tr, there exists a corresponding measurement m = XT{t), where 
'hr = {nT ; n G Z}, and T is a positive real number called the sampling interval. The 
signal domain is dehned as the set If all such measurements are members of the 
extended real numbers, then XTif) is said to be a real-valued (or real) discrete-domain 
signal. If the measurements are members of the complex numbers, then the signal is 
said to be a complex-valued (or complex) discrete-domain signal. When t represents 
time, a discrete-domain signal may be called a discrete-time signal. 

This dehnition of a discrete-domain signal is similar to that of an analog 
signal except that the signal is dehned on a countable set Tr. An important obser¬ 
vation is that a discrete-domain signal is equivalent to a sequence and an associated 
sampling interval T or its reciprocal F = 1/T, 

XT{t) = {x[n],T} = {x[n],F} for n G Z. (2.4) 

The quantity F is called the sampling rate (in samples/sec or Hz) and in discussing 
discrete-domain signals, it is common to refer to the sequence and its sampling rate. 
For example, we may use the expression “a;[n] at a rate of 20 kHz” to describe a 
discrete-domain signal, which has a sampling interval of T = 0.05 msec. 
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It is also common not to mention the sampling rate if the sampling 
rate is common throughout a system (single-rate system). On the other hand, when 
dealing with a multirate system, it is common to use different letters, such as n and m, 
to designate sequences, for example, x[n] and y[m], to indicate that these sequences 
represent discrete-domain signals with different sampling rates. 

Figure 2.1 illustrates a discrete-domain signal. Note that the signal is 
defined only on the points t = nT and is undefined everywhere else. Note, also, that 
while a discrete-domain signal may be derived by sampling an analog signal, this is not 
always the case. Any sequence, regardless of how it is computed (say in MATLAB or 
on an ASIC chip) when combined with a sampling interval, dehnes a discrete-domain 
signal. The corresponding analog signal need not exist unless (as in the output of a 
digital signal processing chain) some special action is taken to construct it. 


xfit) 



Figure 2.1. Graphical representation of a discrete-domain signal XT{t) with sampling 
interval T = 0.05. Note that the signal is defined only at f = nT]n G Z. 

b. Random Signals and Sequences 

In statistical signal processing, a probabilistic model is necessary for 
signals. This model is embedded in the concept of a random signal or a stochastic 
signal. A real random signal or {real stochastic signal) is defined as follows. 
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Definition 4. A real random signal, denoted by or when it is clear from 

context X{t), is a set of ordered random variables (representing measurements) such 
that for every f G M, there exists a corresponding random variable X{t). 

Note that when the context is clear, a random signal may be designated by a lower 
case variable, i.e., x{t), d{t), etc. 

Since a random variable is a mapping from some sample space to the 
real line, the definition for a complex random signal requires special caution. The 
following definition is therefore provided. 

Definition 5. A complex random signal or [complex stochastic signal), denoted by 
{Z[t)}, is defined by Z[t) = X[t)+jY[t), where X[t) and Y[t) are real random analog 
signals defined on a common domain. In other words, for every t G M, there exists a 
pair of corresponding random variables X[t) and Y[t) such that Z[t) = X[t) + jY[t). 
Again, we may use Z[t) instead of {Z[t)} when the meaning is clear from context. 

Random sequences and random discrete-domain signals can be defined 
in a similar manner. 

Definition 6. A real random seguence or [real stochastic seguence), denoted by 
{X[n]}, is a countable set of ordered random variables (representing measurements) 
such that for every n G Z, there exists a corresponding random variable X[n]. A 
complex random seguence can be defined in a manner similar to that of a complex 
random signal. 

Note that when the context is clear, a random sequence may be designated by a lower 
case variable, i.e., x[n], d[n], etc. 

Definition 7. A random discrete-domain signal, denoted by {Xr(t)}, or when it is 
clear from context XT[t), is a set of ordered random variables (representing mea¬ 
surements) such that for every t G there exists a corresponding random variable 
XT[t), where = {nT;n G Z}, and T is the sampling interval. 
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A random discrete-domain signal is sometimes also referred to as a time series; how¬ 
ever, the use of that term in the literature is not always consistent. 

c. Multi-channel Signals and Sequences 
In signal processing, it is often the case that a system may contain 
signals or sequences that are derived from multiple sources or multiple sensors. In 
order to represent such signals and sequences, multi-channel signals and sequences 
are dehned. A multi-channel signal is a set of (single-channel) signals that share a 
common domain and is represented by a vector 

Xi{t) 

X2 

XN(t)^ 

whose components xi(t),X 2 (t),... ,XN(t) are (analog or discrete-domain) signals as 
dehned earlier. The signals may be real or complex, deterministic or random. By 
convention, bold face and vector notation are used to represent such signals as in 


x(t) = 


x(t) 


cos cut 
— sin cut 


or m 


X(() = 


Acos{ut <h) 

—Asin(c<;t -|- <h) 

where X{t) represents a random signal dehned by random variables A and <h. 
A multi-channel sequence 


x[n 


xi[n] 

X2[n] 


XN[n 


17 



is represented by a vector whose components xi[n], X 2 [n],... ,XN[n] are sequences as 
defined earlier. Again, all of the terms describing an individual sequence (e.g., real, 
complex, etc.) can be applied to a multi-channel sequence. 

d. Two-dimensional Signals and Sequences 

Since two-dimensional signals and sequences are at the heart of image 
processing, it is helpful to characterize the 2-D counterparts to the familiar one¬ 
dimensional signals and sequences already presented. A two-dimensional (2-D) analog 
signal is defined as follows. 

Definition 8. A two-dimensional (2-D) analog signal, denoted by {x(ti, t 2 )}, or when 
it is clear from context x{ti, t 2 ), is a set of ordered measurements such that for every 
pair {ti,t 2 ) G there exists a corresponding measurement m = x(ti,t 2 ). Two- 
dimensional signals can be real or complex, deterministic or random. It is sometimes 
convenient to represent a 2-D signal with a bold face argument t = (ti, t 2 ) G Thus, 
the 2-D signal would be denoted by or a;(t) when clear from the context. 

Although a seguence seems to imply an ordered set of terms in one 
dimension, it is common in signal processing to extend the meaning to apply to 
signal defined on a two-dimensional domain. A two-dimensional seguenee and two- 
dimensional discrete-domain signal are thus defined as follows. 

Definition 9. A two-dimensional seguenee, denoted by {x[ni,n 2 ]}, or when it is 
clear from context x[ni,n 2 \, is a set of ordered measurements such that for every pair 
{ni,n 2 ) G I?, there exists a corresponding measurement m = x[ni,n-^. 2-D sequences 
can be real or complex, deterministic or random; they may also be represented as 
{a;[n]} or a;[n], where the boldface argument denotes the ordered pair (ni,n 2 ) G I?. 

Definition 10. A two-dimensional diserete-domain signal, denoted by {xTiT 2 (ti,t 2 )} 
or XTiT 2 (ti, h), is a set of ordered measurements such that for every pair (ti, t 2 ) in the 
domain '^TiT 2 = x ^T 2 i where is as defined earlier, there exists a corresponding 
measurement m = XTiT 2 (ti,t 2 ), and Ti and T 2 are the associated sampling intervals. 
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For convenience in notation, we may use a;T(t) and 'Ft to denote the 2-D signal 
and its domain, where T represents the ordered pair (Ti,T 2 ) of sampling intervals. 
Again, note that a two-dimensional discrete-domain signal can be real or complex, 
deterministic or random. 

The image projected on the him plane of a camera is an example of 
a 2-D analog signal. If him is thought of as a continuous medium, then the image 
captured on the him is also a representation of a 2-D analog signal. If the image is 
projected onto a sensor array as in a digital camera, then the resulting sampled image 
is represented by a 2-D discrete-domain signal. 

Signals can be both multi-dimensional and multi-channel. A common 
example is a color image where the domain is two-dimensional (horizontal and vertical 
spatial variables), and there are 3 channels corresponding to the three components of 
a color space, such as RGB (red, green, blue), CMY (cyan, magenta, yellow) or HSI 
(hue, saturation, intensity). 

Two-dimensional random signals and sequences are similar to their cor¬ 
responding deterministic representations except that the measurements are repre¬ 
sented by random variables. 
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e. Summary of Notation and Convention 

A summary of the various signal representations is provided in Ta¬ 
ble 2.1. 


Representation 

Name 

x{t) 

Deterministic analog signal, analog aignal 

x[n] 

Deterministic sequence 

XT{t),x[n]T 

Deterministic discrete-domain signal with sampling in¬ 
terval T, Discrete-domain signal 

x{ti,t- 2 ),x{t) 

Two-dimensional deterministic analog signal, 2-D analog 
signal 

a;[ni,n2],a;[n] 

Two-dimensional deterministic sequence, 2-D determin¬ 
istic sequence 

XTiT2{iii h), a^T(t) 

Two-dimensional deterministic discrete-domain signal 
with sampling intervals Ti and T 2 , 2-D discrete-domain 
signal 

X{t) 

Random analog signal 

X[n] 

Random sequence 


Table 2.1. Summary of signal representations. 


B. CONCEPTS IN LINEAR ALGEBRA 

1. Random Vectors 

Often, it is necessary to process some finite number of samples of a random 
sequence. Such a hnite-length sequence can be conveniently represented by a random 
vector [Ref. 5]. This provides for compact notation and formulation and solution of 
problems in a linear algebra sense. A random sequence X[n] restricted to some 
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interval 0 < n < — 1 can be represented by an iV-component random vector x as 

shown in Figure 2.2 and written as 


A'|0] 



Figure 2.2. Graphical representation of a finite-length random sequence as a random 
vector. 


2. Kronecker Products 

The Kronecker product, also known as the direct product or tensor product, 
has its origins in group theory [Ref. 4] and has important applications in a number of 
technical disciplines. In this study, the Kronecker product is used to develop matrix 
representations of various multirate operations. 


Definition 11. Let A be an m x n matrix (with entries aij) and let B be an r x s 
matrix. Then the Kronecker product of A and B is the mr x ns block matrix 

A 


A G) B = 


f aiiB 

ai2B . 

dinB 

a2iB 

n22B 



®m2B 



( 2 , 6 ) 


Equation (2.6) is also called a right Kronecker product as opposed to the 
definition A (g)' B = B (g) A, which is called a left Kronecker product. Since there is 
no need to use both, we will stick with the more common dehnition (2.6). 
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A summary of some important properties of the Kronecker product is provided 


in Table 2.2. 


A ® (q:B) = a(A ® B) 

(A + B)®C = A®C + B(8)C 
A ® (B 0 C) = (A ® B) ® C 
(A ® B)^ = A^ ® B^ 

(A ® B)(C ® D) = AC ® BD 

(A® B)-i = A-i ®B-i 


Table 2.2. Some Kronecker product properties and rules, (After [Ref. 4]). 


3. Reversal of Matrices and Vectors 

In signal processing, it is a common requirement to view signals as evolving 
either forward or backward in time. A well-known example is the convolution opera¬ 
tion, where the linear combination of terms involves a time-reversed version of either 
the input signal or the system impulse response. Since, in discrete-time signal pro¬ 
cessing, signals are often represented by vectors, it is useful to dehne the operation 
of reversal for vectors and matrices. 

The reversal of a vector x is the vector with its elements in reverse order. 
Given the vector 



Xi 


Xn 

X = 

X-2 

, its reversal is x = 

Xn-1 


Xn 


Xi 


Note that the notation for the reversal is x, and it is used just like notation for the 
transposition of a vector or matrix. 
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The reversal of a matrix A is the matrix with its column and row elements in 


reverse order. Given the matrix A G M 

MxN 




ail 

0^12 


OlN 

A = 

^21 

0-22 


02N 


ttMl 

0M2 


omn 

its reversal A G is given by 





O-MN 


0M2 

omi 

A = 

0-2N 


022 

021 


aiN 


Ol2 

Oil 


Note that the reversal of a vector or matrix can be formed by the product of a 
conformable counter identity and the vector or matrix itself. 

Some common properties of the reversal operator are included in Table 2.3. 
In particular, the reversal of matrix and Kronecker products (see Section II.B.2) are 
products of the reversals, and the operation of reversal commutes with inversion, 
conjugation and transposition. 



Qnantity 

Reversal 

Matrix prodnct 

AB 

AB 

Matrix inverse 

A-i 

(A)-i 

Matrix conjngate 

A* 

(A)* 

Matrix transpose 

A^ 

(A)^ 

Kronecker prodnct 

A (g) B 

A (g) B 


Table 2.3. Some properties of the reversal operator, (After [Ref. 5]). 
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4. Probenius Inner Product 

In the development of approaches to two-dimensional signal reconstruction, it 
is convenient to express the related linear estimates in terms of the Frobenius inner 
product. 

Definition 12. For any A, B G with elements aij,bij, the Frobenius inner 

product of the matrices is dehned as 

m n 

(A, B) = tr(AB^) (2,9) 

1=1 j=l 

C. MOMENT ANALYSIS OF RANDOM PROCESSES 

Generally, a complete statistical model is unavailable when analyzing systems 
of random processes. Either the required joint density functions are unavailable, or 
they are too complex to be of utility. If the random processes under consideration are 
Gaussian, then the system can be fully specihed by only its hrst two moments [Ref. 5]. 
Even if the processes are not Gaussian, second moment analysis is often adequate in 
characterizing the statistical relationships between signals in such systems and forms 
the basis for any additional analyses. This section introduces the required dehnitions 
and relevant properties associated with second moment analysis [Ref. 5]. 

1. Definitions and Properties 

Given the random process A[n], the first moment or mean of the random 
process is dehned by 

mx[n] = £{X[n]}, (2.10) 

where £{■} denotes expectation. 

The correlation between any two samples of the random process X[ni] and 
X[no] is described by the correlation function or autocorrelation function, which is 
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defined by 


Rx[ni,no] = £{X[ni]X*[no]}. (2.11) 

In certain applications, and extensively in this work, it is convenient to define 
a time-dependent correlation function as 

Rx[n-, 1] = 8,{X[n]X*[n - /]}, (2.12) 

and the various dehnitions and relationships introduced in this section will be based 
on this “time-dependent” representation. 

The covariance between any two samples of the random process X[n] and 
X[n — 1] is described by the time-dependent covariance function, which is dehned by 

Cx[n', 1] = £{(X[n] — mx[n]){X[n — /] — mx[n — /])*}• (2-13) 

The relationship between the correlation function and the covariance function is 

Rx[n-, 1] = Cx[n-, 1] + mx[n]m*x[n - /], (2.14) 

hence when X[n] is a zero-mean random process, 

Rx[n-, 1] = Cx[n-, /]. 

If we consider two random processes, X[n] and Y[n], the correlation between 
any two samples of the random processes is described by the time-dependent cross¬ 
correlation function, which is dehned by 

RxY[n-,l] = E{X[n]Y*[n-l]}. (2.15) 

An expression can be written for the time-dependent cross-covariance function as 

Cxvin; 1] = £{(X[n] — mx[n]){Y[n — /] — my[n — /])*}. (2.16) 

The relationship between the cross-correlation function and the cross-covariance func¬ 
tion is 

RxY[n-, 1] = CxY[n-, 1] + mx[n]my[n - /], (2.17) 
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hence when X[n] and Y[n] are zero-mean random processes, 

RxY[n'-,l] = CxY[n'-,l]- 

Two random processes are called orthogonal if RxyW^I] = 0 and uncorrelated if 
/] = 0 . 

2. Stationarity of Random Processes 

Recall that a random process is wide-sense stationary (WSS) if 

1. the mean of the random process is a constant, mx[n\ = mx, and 

2. the correlation function is a function only of the spacing between samples, i.e., 
Rx [n; 1] = Rx [1] ■ 

and that two random processes are jointly wide-sense stationary (JWSS) if 

1. they are each WSS, and 

2. their cross-correlation function is a function only of the spacing between sam¬ 
ples, i.e., RxY[n;l] = Rxy[1]- 

Under the assumptions of WSS and JWSS, the mean, correlation and covari¬ 
ance functions are summarized in Table 2.4. 

3. Matrix Representations of Moments 

Using the vector representation (2.5) for a random signal, a number of impor¬ 
tant concepts and properties can be defined. The first moment or mean of a random 
vector is defined by 

£{X[0]} j r mx[0] 

£{X[1]} mx[l] 

mx = c|X| = = _ 

£{X[iV-l]} mx[N-l] 
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Mean Function 

mx = £{X[n]} 

(2.18) 

(Auto)correlation Function 

Rx[l] = e,{X[n]X*[n - 1]} 

(2.19) 

Covariance Function 

Cx[l] = £{(X[n] - mx){X[n - 1] - mx)*} 

(2.20) 

Interrelation 

Rxi^] = Cx[l] + \mx\'^ 

(2.21) 

Cross-correlation Function 

RxY[l] = ^{X[n]Y*[n-l]} 

(2.22) 

Cross-covariance Function 

Cxy[1] = £{(X[n] - mx)(Y[n - 1] - my)*} 

(2.23) 

Interrelation 

Rxy[1] = Cxy[1] + mx^riy 

(2.24) 


Table 2.4. Summary of definitions and relationships for stationary random processes, 
(After [Ref. 5]). 


which is completely specihed by the associated mean function mx[n] in (2.10). If the 
random process is WSS, then the mean function is independent of the sample index 
and mx is dehned by a vector of constants 


mx 


mx 


mx 


(2.26) 


mx 


The correlation matrix represents the complete set of second moments for the 
random vector and is dehned by 


Rx = £{XX*^}. 


(2.27) 
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The correlation matrix thus has the explicit form 


Rx = 


£{|X[0]p} 

£{X[1]X*[0]} 


£{X[0]X*[1]} 

£{|X[l]p} 


£{X[0]X*[iV- 1]} 
£{X[l]X*[iV- 1]} 


(2.28) 


£{|X[iV- 1]X*[0]} 

Rx[0;0] 

-Rx[l; 1] 


£{X[iV- 1]X*[1]} 

-Rx[0;—1] 
i?x[l;0] 


. £{|X[iV-l]p} 

i?x[0; —N + 1] 
Rx[l-,-N] 


(2.29) 


Rx[iV-l;iV-1] Rx[N-1;N-2] 


Rx[N — 1; 0] 


which is completely specihed by the associated correlation function Rx[n', 1] in (2.12). 
If the random process is WSS, then the correlation is a function of only the sample 
spacing and has the form of a Toeplitz matrix; 


Rx [0] 

Rx[-l] 

Rx[-2] ... 

Rx[-N+1] 


Rx[l] 

Rx [0] 

Rx[-1] 

Rx[-N] 


Rx[2] 

Rx[l] 



(2.30) 

Rx[N-l] 

Rx[N-2] 

... Rx[l] 

i?x[0] 



This matrix is completely specihed by the associated correlation function Rx[l] in 
(2.19). 

The cross-correlation matrix represents the complete set of second moments 
between two random vectors X G and Y G and is dehned by 


Rxy = £{XY'’’}, 


(2.31) 




and the associated correlation matrix has the form 


-Rxr [0; 0] Rxy[0', —1] • • • -Rxr [0; —M + 1] 

-Rxy[l;l] -Rxy[l;0] ... Rxy[^',—M] 

Rxy = , 

Rxy[N - 1 ] N - 1 ] Rxy[N - 1 ] N - 2] ... Rxy[N - 1 ] N - M] 

(2.32) 

which is completely specified by the associated cross-correlation fnnction RxY[n; 1] in 
(2.15). In general, Rxy is not a sqnare matrix (unless N = M). If the associated 
random processes are JWSS, then the cross-correlation is a function of only the sample 
spacing 

Rxy[0] Rxy[-1] ... Rxy[-M+1] 

Rxy[1] Rxy[0] Rxy[-M] 

Rxy = Rxy[‘^] Rxy[1] '. ... > (2.33) 

_Rxy[N-1] Rxy[N-2] ... Rxy[N-M]_ 
which is completely specihed by the associated correlation function Rx[l] in (2.22). 
In general, such matrices will exhibit Toeplitz structure but will not be Hermitian 
symmetric [Ref. 5]. Similar expressions and statements can be made concerning 
the cross-covariance matrix and function. The essential dehnitions, properties, and 
relations for the quantities discussed in this section are listed in Table 2.5. 

4. Reversal of First and Second Moment Quantities 

Since the operations of expectation and reversal commute, we have the follow¬ 
ing relations for the first and second moment quantities 

= £{X} = mx, (2.34) 

and 

Rx = £{XX*^} = Rx (Cx = Cx). (2.35) 
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Mean 

mx = £{X} 

(Auto) correlation 

Rx = £{XX*^} 

Covariance 

Cx = £{(X-mx)(X-mx)*^} 

Interrelation 

Rx = Cx + mxmx*^ 

Cross-correlation 

Rxy = £{XY*^} 

Cross-covariance 

Cxy = £{(X - mx)(Y - itiy)*^} 

Interrelation 

Rxy = Cxy + mxniY*^ 

Symmetry 

ivx — ivx , '.-^X — '.-"X 

Relation of Rxy and Cxy 

ivxY — iVYX , '.-"XY — '.-"YX 


Table 2.5. Summary of useful definitions and relationships for random processes, 
(After [Ref. 5]). 

Further, if Rx (Cx) is a Toeplitz correlation (covariance) matrix corresponding to a 
WSS random process, it follows that 

R;, = Rx. (2.36) 


D. NUMBER THEORY 

Number theory, “... the branch of mathematics concerned with the study of 
the properties of the integers [Ref. 55],” is a natural framework for the analysis of 
discrete-time systems, where the independent variables, by dehnition, are integers. 
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In particular, since in this analysis of multirate systems, notions of divisibility, factor¬ 
ization and congruence are integral, the ensuing discussion is provided to introduce 
and dehne these and related concepts [Ref. 55, 56, 57, 58]. 

1. Division Algorithm Theorem 

The elementary operation of division forms the basis of much of what is to 
follow and is expressed by the division algorithm theorem. 

Theorem 1. Let a and b be integers with a > 0. Then there exists unique integers 
q and r satisfying 

b = qa + r, 0 < r < a, (2.37) 

where q is called the quotient and r is called the remainder. 

The proof of this can be found in many texts, e.g., [Ref. 55, 56, 57]. 

Example 1. A specifie example demonstrating the division algorithm theorem is pro¬ 
vided. Given integers a = 3 and b = 22, we find unique integers q = 7 and r = 1 that 
satisfy (2.37). The quotient is q = 7; the remainder is r = 1. ■ 

2. Divisibility 

Definition 13. Let a and b be integers. Then a divides b, written a\b, if and only if 
there is some integer c such that b = ca. When this condition is met, the following 
are equivalent statements: (i) a is a factor of b, (ii) b is divisible by a, and (iii) 5 is a 
multiple of a. If a does not divide b, we write a\ b. 

Example 2. This example illustrates the coneept of divisibility for a number of integer 
pairs. 

3|12,7|21,9|108,12|144; 

4 {5, 7 {8, 8 {7, 3 {22. ■ 
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a. Greatest Common Divisor 

Definition 14. Let a and b be integers. The integer d is called the greatest eommon 
divisor of a and h, denoted by gcd(a, b), if and only if 

1. d > 0, 

2. d\a and d\b, 

3. whenever e|a and e\b, we have e\d. 

The integers a and b are said to be relatively prime if gcd(a, b) = 1. 

Example 3. A few examples demonstrating the greatest eommon divisor: 

If a = 3 and 6 = 4, then d = gcd(3,4) = 1 (3 and 4 are relatively prime), 

If a = 12 and b = 15, then d = gcd(12,15) = 3, 

If a = 25 and b = 55, then d = gcd(25, 55) = 5. ■ 

b. Least Common Multiple 

Definition 15. Let a and b be positive integers. The integer m is called the least 
eommon multiple of a and b, denoted by lcni(a, b), if and only if 

1. m > 0, 

2. a\m and b\m, and 

3. if n is such that a|n and b\n, then m\n. 

The least common multiple can be expressed as 

nh 

lcm(a,6) = (2.38) 

gcd(a, b) 

Example 4. A few examples demonstrating the least eommon multiple: 

If a = 3 and 6 = 4, then m = lcm(3,4) = 12, 

If a = 12 and 6 = 15, then m = lcm(12,15) = 60, 

If a = 25 and 6 = 55, then m = lcm(25, 55) = 275. ■ 
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Also note that the least common multiple is associative and therefore, 

lcm(a, b, c) = lcm(lcm(a, b),c) = lcm(a, lcm(6, c). (2.39) 

3. Greatest Integer Function 

The greatest integer funetion, often called the floor funetion, is dehned as 
follows. 

Definition 16. For any a; G M, the greatest integer function evaluated at x returns 
the largest integer less than or equal to x. This is sometimes referred to as the integral 
part of X. The function will be denoted as [xj. 

Example 5. The following examples illustrate this definition, 

L2.7J = 2, 

L0.9J = 0, 

[-0.3J = -1. ■ 

Note that the floor function satishes the following identity 

[x + k\ = [xj + k, for /c G Z. (2.40) 

4. Congruence 

If a is hxed in (2.37), then there are an inhnite number of choices of b for which 
the remainder r is the same. In this context, a is called the modulus, the choices of b 
are said to be congruent modulo a, and the remainder is called the common residue 
modulo a or simply the common residue [Ref. 58]. This concept of congruence is 
formalized with the following dehnitions. 

Definition 17. Let n be a positive integer. The integers x and y are “congruent 
modulo n” or “x is congruent to y modulo n”, denoted x = y (mod n), provided that 
X — y is, divisible by n. If x and y are not congruent modulo n or a; is not congruent 
to y modulo n, we write x ^ y (mod n). 
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Example 6. We demonstrate the coneept of congruenee with a few examples. 

8 = 5 (mod 3), 

14 = 2 (mod 12), 

49 = 42 (mod 7). ■ 

Example 7. In the following example, n = 2, and there are two sets of integers that 
are congruent modulo 2, the even integers and the odd integers. 

, —4, —2, 0, 2,4,.. .} are congruent to 0 (mod 2), 

, —3, —1,1, 3, 5,.. .} are congruent to 1 (mod 2). ■ 

Example 8. In this example, n = 3, and there are three sets of integers that are 
congruent modulo 2. 

, —6, —3, 0, 3, 6, ...} are congruent to 0 (mod 3), 

, —5, —2,1,4, 7, ...} are congruent to 1 (mod 3), 

, —4, —1, 2, 5, 8, ...} are congruent to 2 (mod 3). ■ 

Definition 18. If x = y (mod n), then y is called a residue of x modulo n. Further¬ 
more, if 0 < 1 / < n, then y is called the common residue of x modulo n, or simply the 
common residue. 

Example 9. Referring to Example 6, we point out the associated residues. 

5 is a residue of 8 modulo 3, 

2 is the common residue of 14 modulo 12, and 
42 is a residue of 4Q modulo 7. ■ 

Definition 19. The set of integers A„ = {0,1,... ,n — 1} is called the set of “least 
positive residues modulo n” 

At times, it is necessary to extract the common residue [Ref. 58]. This oper¬ 
ation is denoted by (■)„ and is defined as 

X 

y = {x)n = X- - n, (2.41) 

n 

where y is the common residue of x modulo n, and [-J is the floor operation. 
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Example 10. A few examples of extracting the common residues of x modulo n. 


y = (22)3 
y = ( 14)4 


22- 


14- 


3 


3 = 22-21 = 1, 


4 = 14- 12 = 2. 


E. CHAPTER SUMMARY 

This chapter introduces various fundamental signal processing and mathemati¬ 
cal concepts required for theoretic and application-related developments in subsequent 
chapters. Further, for the purposes of consistency, accuracy and ease of understand¬ 
ing, conventions and notation are also established. 

The taxonomy of signals and sequences, their various definitions, and associ¬ 
ated notations are presented. Of particular relevance is the discussion on discrete- 
domain signals and their sequence representation, which form the most basic con¬ 
stituent of any multirate system (Chapter III). 

Many concepts from linear algebra are recalled, including the concept of a 
random vector and the reversal of a vector or matrix. Further, the linear algebraic 
concept of the Kronecker product is discussed, which is useful in the matrix represen¬ 
tation of various multirate operations in Chapter III and the multirate Wiener-Hopf 
equations in Chapter IV. Finally, the Frobenius inner product is introduced, which 
provides a compact representation of the two-dimensional linear estimate required for 
image reconstruction (Chapter V). 

In the analysis of random processes, the second-moment properties are fre¬ 
quently used. Since they are essential to the development of optimal estimation 
theory, the analysis and various definitions and relationships are reviewed in this 
chapter. 


35 



Finally, several topics in number theory are presented, which have great utility 
in developing the theory of multirate systems and characterizing the relationships 
between constituent signals and the related multirate system (Chapter III). 
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III. MULTIRATE SYSTEMS: CONCEPTS AND 

THEORY 

In this chapter, we develop the theory of multirate systems, which establishes 
the fundamental relationships in a multirate system, and culminates in a systematic 
framework for their analysis. These results lead to representation of the various 
signals in a multirate system on a common domain, system and impulse response 
formulations at both the signal- and system-level, linear algebraic representation of 
multirate operations, and ultimately, as presented in Chapter IV, development of 
multirate signal estimation theory. 

A. INTRODUCTION 

In many digital signal processing (DSP) applications, the systems involved 
must accommodate discrete-domain signals that are not all at the same sampling 
rate. For instance, consider a system in which the signals at the source and desti¬ 
nation have different sampling rate requirements. An example of this occurs when 
recording music from a compact disc (CD) system at 44.1 kHz to a digital audio tape 
(DAT) system at 48 kHz. Another application might involve systems that incorporate 
several signals collected at different sampling rates. Sensor networks, many military 
weapon and surveillance systems, and various controllers process data from various 
sensors, where the information from each sensor might be collected at a different rate. 
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Further, a system may be at a rate that is inefficient, and sampling rate conversion 
may be required to reduce the rate of the system because “oversampling” is wasteful 
in terms of processing, storage and bandwidth. 

B. MULTIRATE SYSTEMS 

The various ideas described in this chapter follow [Ref. 3, 59]; however, many 
important extensions are made to align results with the theory of multirate systems as 
developed here. A multirate system will be defined as any system involving discrete- 
domain signals at different rates. Recall from Chapter II that we will use sequence 
notation (i.e., x[n]) and different index values {n, m, etc.) to denote discrete-domain 
signals at different rates. Figure 3.1 depicts a notional multirate system where the 
input, output and internal signals are at different rates. A specihc example of a mul- 



Figure 3.1. Notional multirate system where input, output, and internal signals are 
at different rates. (From [Ref. 3]). 

tirate system is the subband coder illustrated in Figure 3.2. The signals x[n] and y[n\ 
at the input and output of the system are at the original sampling rate while some of 
the internal signals {yi[mi\ and 2 / 2 [^^ 2 ]) are at lower rates produced through hltering 
and decimation. 
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G2{z) 


y^im^] z^[n] 

Figure 3.2. Simple subband coding system. 


1. Intrinsic and Derived Rate 

The notion of rate was introduced in Chapter II and is part of the description 
of any discrete-domain signal. The rate associated with a particular signal may be 
a result of sampling an analog signal or a result of operations on sequences in the 
system. These issues are discussed below. 

a. Intrinsic Rate 

A discrete-domain signal may be derived from an analog signal by pe¬ 
riodic or uniform sampling described by 

x[n] = x{nT^) = x{t)\t=nT^ n E Z. (3.1) 

Here, x[n\ is the discrete-domain sequence obtained by sampling the analog signal 
x{t) every seconds. This concept is depicted in Figure 3.3. 


x{t) 

— V— 

x[«] = x{nTJ 


T 

X 



Figure 3.3. An analog signal sampled with a sampling interval of T^. 
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The sampling interval and its reciprocal, the sampling rate are related by 

F. = (3.2) 

In this context, we say x[n] is at a rate F^. The rate associated with the sequence, 
therefore, is the rate at which its underlying analog signal was sampled and is referred 
to as its intrinsic rate. 

b. Derived Rate 

The process of sampling rate conversion provides another context in 
considering the notion of rate or sampling rate in multirate systems. The two basic 
operations in sampling rate conversion are downsampling and upsampling (with ap¬ 
propriate hltering). These operations are depicted by the blocks shown in Figure 3.4, 
and they are mathematically represented by 

y[n]=x[Mn], (3.3) 


where n is an integer, in the case of downsampling, and 

[x [fl , if n\L ; 

y[n]=l (3.4) 

I 0, otherwise, 

in the case of upsampling. Figures 3.5 and 3.6 graphically depict the downsampling 
and upsampling operation, respectively, for M = L = 2. 



Figure 3.4. Basic operations in multirate signal processing, downsampling and up- 
sampling. 
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Figure 3.5. An example of the downsampling operation (3.3), M = 2. 



Figure 3.6. An example of the upsampling operation (3.4), L = 2. 
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Note that both operations are performed exclusively in the digital domain. The 
resulting signals have no intrinsic rate, but the rate is derived from the rate of the 
input signal. For downsampling, the output rate Fy is given by 


F - IT 
" M’ 


while for upsampling the rate is given by 


(3.5) 


Fy — LFx. (3.6) 

The parameter M in downsampling is called the decimation factor while the parame¬ 
ter L may be called the upsampling factor. Thus, downsampling results in a reduction 
of sampling rate by a factor of M, and upsampling results in an increase in sampling 
rate by a factor of L. 

It will be seen later that other operations more general than downsam¬ 
pling and upsampling can result in rate changes. These more general operations will 
be represented by linear periodically-varying hlters (see Section III.D.2). The out¬ 
puts of these hlters have no intrinsic rate but have a derived rate associated with the 
operation that is performed. 


C. CHARACTERIZATION OF MULTIRATE SYSTEMS 

In the discussion of multirate system concepts and associated theory, it is 
necessary to further develop terminology, characterize such systems and develop a 
conceptual framework by which further analysis and extension can be based. In this 
section, the concepts and terms are introduced. 

1. System Rate 

Consider a multirate system with just two signals at different sampling rates, 
Xi and X 2 - Although it is not strictly necessary, the discussion can be more easily 
motivated if it is assumed that each signal is derived by sampling an underlying analog 
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signal as shown in Fignre 3.7. It will be assnmed that the sampling rates Fi and F 2 
are integer-valued. While the treatment could be generalized to the case where the 
rates are rational numbers, the assumption of integer-values simplihes the discussion 
and is quite realistic for practical systems. The corresponding discrete-domain signals 


Vi(0 

— V— 






^2(0 




— 

F2 



Figure 3.7. Two signals sampled at different sampling rates. 


xti and xt 2 at the output of the samplers are defined at points on their respective 
domains 



= {nTi ;n G Z}, 

(3.7) 

and 

4/^2 = {nT 2 ;n e Z}, 

(3.8) 

where Ti = — and T 2 = —. The discrete-domain signals are represented in Figure 
Fi F 2 

3.8 as sequences with different index values x[n]\ and x[n 2 \ indicating the different 

sampling rates. 

Note that there is some common domain 



Tjt = \nT ; n G Z}, 

(3.9) 


with some maximum sampling interval T in which the samples of both xi and X 2 can 
be represented. In other words, Tti C and C Tyr. 

The sampling interval T in (3.9) will be called the system sampling interval or clock 
interval. We can state the following theorem. 
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Figure 3.8. Two signals sampled at different integer-valned sampling rates. A periodic 
correspondence between indices can be observed (as indicated by the dashed lines). 
The system grid is represented by the line segment at the bottom of the hgnre and 
is derived from the set of hidden and observed samples of the associated underlying 
analog signals. Open circles represent “hidden” samples. 


Theorem 2. The system sampling interval is given by T 
integer, and 

F = lcm(Fi, F 2 ). 

Here, F is called the system rate or the fundamental rate. 

Proof. From the dehnition of the least common mnltiple (Definition 15): 

1. By dehnition, F must be a positive integer; 

2. Tti C % ^ Fi\F, ^t2 C % ^ F 2 IF; 

3. Since T is the maximnm sampling interval in which samples of both xi and 
X 2 can be represented in this implies that F = = is the related minimum 
sampling rate, and the third condition of Dehnition 15 is met. 


= =, where F is a positive 


(3.10) 


□ 
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We also define a system grid as the representation of the set on the real line, 
as seen in Figure 3.8. Note that the samples of a signal corresponding to the system 
grid can be viewed as a set of hidden and observed samples of the underlying analog 
signal. For example, for a given signal xi, the observed samples are the samples of 
xi(t) that correspond to the set d'ri dehned in (3.7) while the hidden samples of xi 
are the samples of xi(t) that correspond to the set 

It is frequently useful to represent all signals at the “system level” dehned 
by T and F. The system level is referred to as the “fundamental layer ” in [Ref. 
59]. Sequences associated with the system level will be denoted by a symbol with an 
overbar, as in xi [n] and X 2 [n] and a common index n. This point is developed further 
in Section III.C.5. 

For a multirate system comprised of M signals, the dehnition (3.10) can be 
extended (see (2.39)) as 

F = \cm{Fi,F2,.. . ,Fm), (3.11) 

with T = =. 

F 

2. Decimation Factor 

Recall, from Section II.D.2(b), that the least common multiple m is a number, 
which is a multiple of both associated integers a and b, therefore, a|m and b\m. 
Further, from Dehnition 13, the condition a|m implies that m = cia, and likewise 
b\m implies that m = C 26 , where ci and C 2 are constants. 

We can apply these results to multirate systems. Since F is the least common 
multiple of Fi and F 2 , it follows that there exists integer constants Ki and K 2 such 
that 

F = KiFi and F = K 2 F 2 . (3.12) 

^The “difference” A — B of two sets A and B, where B C T, is defined as A — B = {x € A\x ^ B}. 
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Rearranging Equation (3.12) yields an expression for these constants as 

R'l = 4" and K 2 = 4-- (3.13) 

J^i r2 

Further, for a multirate system comprised of M signals, (3.13) can be extended as 


F 

Ki = — where i = 1 ,..., M. 
Fi 

Notice that if (3.14) is rearranged to the form of (3.5): 

F 


(3.14) 


Fi = 


Ki' 


we explicity see that the system sampling rate F is reduced by a factor of K^. There¬ 
fore, Ki is dehned as the system decimation factor, or simply the decimation factor, 
for the signal. 

As a consequence, the following observation can be made and extended to 
any number of signals. If the underlying analog signals xi and X 2 are sampled at 
the system rate F and then are decimated by their respective decimation factors Ki 
and K 2 , the resultant discrete-domain signals are the signals obtained by sampling 
xi and X 2 at the original rates Fi and F 2 , respectively (see Figure 3.9). In terms of 
the hidden and observed samples of xi and X 2 , the decimation factors offer a way to 
relate the set of observed samples to the set of all samples (hidden and observed) at 
the system level. 

If (3.14) is expressed in terms of sampling intervals, we have 

(3.15) 

where the decimation factor is the ratio of the duration of the signal’s sample 
interval to the duration of its associated clock interval. Thus, the decimation factor 
Ki can also be viewed as the number of clock intervals or system samples between 
samples of the signal when sampled at T). 


T- 

Ki = =, 

T 


3. System Period 

Again, consider a multirate system comprised of two signals xi and X 2 sampled 
at integer-valued sampling rates Fi and F 2 , respectively. If samples of these signals 
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Figure 3.9. Signals sampled at the system rate and decimated by their respective 
decimation factors yield the original discrete-domain signals. 

align on the related system grid at some time t = to, then their samples realign 
at some later time t = ti = Iq + T and at regular intervals thereafter (see Figure 
3.10). The smallest positive value of T for which the realignment occurs is called 
the system period (not to be confused with the system sampling interval T).Notice 
that T = ti — to and if ti and to are expressed in terms of their discrete-domain 
representations, we can write T = (ni — no)T; thus, we dehne 

N = ni- no, 

which is called the discrete system period. Observe that N represents the number 
of system samples between sample realignments. It relates to the system period (in 
seconds) as 

T = NT. (3.16) 


Now, we dehne Mi to be the number of signal samples per period. Therefore, 
we relate Mi and M 2 to the discrete system period by MiTi = NT = M 2 T 2 . Recalling 
(3.14), we write 

N = MiKi = M 2 K 2 , (3.17) 

and recalling Dehnition (II.G.b.15), we can write 

N = \cm{Ki,K2), (3.18) 
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x,{t) 



Figure 3.10. Two signals sampled at different integer-valued sampling rates. Observe 
the periodic alignment between indices, (After [Ref. 3]). 

an explicit expression for the discrete system period. Also notice from (3.17) that 

N N 

Ml = — and M 2 = —. (3.19) 

Ki K 2 

For a multirate system comprised of M signals, the dehnition of (3.18) can be 
extended (see (2.39)) as 

N = \cm{Ki,K2,...,KM). (3.20) 

4. Maximally-decimated Signal Set 

Consider a discrete-domain signal y[n] at rate Fy. If y[n] is downsampled by 

a factor of M, after M — 1 successive translations, a set of related discrete-domain 

signals results, designated {xo[m],xi[m],... ,XM-i['nT]}. Note that the constituent 

p 

discrete-domain signals are at rate ~ example is shown in Figure 3.11 for 

M = 3. If the discrete-domain signals associated with y[n] are described by 

Xi[m] = y[i + riM] i = 0,1,..., M — 1, (3.21) 
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Figure 3.11. Example of a 3-fold maximally decimated signal set. 


then the resultant signals will be called a maximally-decimated signal set with down- 
sampling factor M or an M-fold maximally-decimated signal set. 

5. Representation of Signals in Multirate Systems 

In the analysis of multirate systems, there is a need to relate signals to a 
common scale. This scale is represented by the system grid and is discussed in terms 
of the domain, where T is the system sampling interval. The sample indices of a 
signal sampled at the system rate correspond to the integer multiples of T. In this 
context, every signal in a multirate system can be represented on the system grid. 

Consider a multirate system, with a system sampling interval T, containing a 
particular discrete-domain signal xt^, among others. This signal can be represented 
by the sequence 

x[mx\ = XT,{mxT^), 

where xt^ is the discrete-domain signal and is the associated sampling interval. 
Note that the sequence and its sampling rate is dehned by the discrete-domain signal 
XT^{t) for t = nixTx and not by the analog signal x(t), which may not be directly 
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observable or may not even exist! We refer to this as a signal-rate or signal level 
representation of x and, in this case, say that x is at its native rate. Now dehne x^ 
to be the associated system-rate or system-level representation of x. This signal can 
be represented by the sequence 

x[n] = x^inT). 

We can relate the two representations by recalling from (3.15) that = TK^ and 
noting that 

xtS^xTx) = xtS^xTK^) = Xj^{mxK^T) =x[K^mx]. 

Therefore, a discrete-domain signal at its native rate can be represented at the system 
rate or system level by 

x[mx] =x[K^mx], (3.22) 

where is the decimation factor of signal x. These concepts are illustrated in 
Figure 3.12 for a multirate system comprised of two signals x and y. The native rate 
corresponds to a signal’s “original” sampling rate, and the system rate corresponds 
to the system rate determined from (3.10). 

6. Summary of Multirate Relationships 

For convenience, the representations and relationships developed in this chap¬ 
ter regarding multirate systems are summarized in the following tables. In Table 
3.1, the various signal representations, their notations and their associated rates are 
indicated. 


x(t) Analog 

x[n] System-level, system rate, F 

x[mx] Signal-level, native rate, 

Table 3.1. Signal representations in multirate systems. 
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Discrete-Time 


x[njJ = Xj. {rnj^) Native-rate 


x[n] = x(Tn^ 

System-rate 

y[n] = y(Tn') 


y[m,]=yr,(m^T^) Native-rate 


Figure 3.12. Signal representations and sampling levels in a multirate system. 


In Table 3.2, the characterization of and relationship between signals in multirate 
systems is summarized. 


Name 

Relationship 


System Rate 

F = lcm(Fi, F 2 , .. 

■, Fm) 

System Sampling Interval 

1 

f _ 


Decimation Factor 

Ai = F = 2 

Fi T 


System Sample Period 

N = \cm{Ki,K2,. 


Samples/Period 



System Period 

T = NT = MiTi 



Table 3.2. Summary of various relationships pertaining to a multirate system (M 
signals). 
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Finally, in Table 3.3, these characterizations and relationships are summarized in 
relation to signal and system level representations. 



Sampling 

Sampling 

Decimation 

Samples 

System 


Rate 

Interval 

Factor 

per Period 

Period 

Signal 

F. 

u 

K. 

M, 

T = 

System 

F 

T 

1 

N 

T = NT 


Table 3.3. Parameters pertaining to a multirate system, (After [Ref. 3]). 


D. MULTIRATE SYSTEM THEORY 

1. Description of Systems 

A signal processing system represents the process for transforming a signal 
into another signal. The concept is illustrated in Figure 3.13(a), where a signal x{t) 
is transformed into a signal y{t). The input and output may be of any of the signal or 
sequence types discussed in Chapter II and need not be of the same type. For instance, 
the input may be a discrete-domain signal and the output may be an analog signal. 
Our primary concern, however, is the case where both the input and the output are 
discrete-domain signals, not necessarily with the same sampling intervals. 

Since discrete-domain signals are represented by sequences together with a 
(known or implied) sampling rate, it is appropriate to consider the properties of 
“systems” that transform sequences. These will be referred to as discrete systems 
as shown in Figure 3.13(b), and the transformation will be represented by 7d- In 
addition to performing a specihed operation on the input sequence, a discrete system 
provides a means of determining the output sampling rate from the input sampling 
rate. When input and output sampling rates are not the same, it is our convention 
to use different letters for the index variable of the sequence. 
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x{t) 


xo 

-► 




r{) 


(a) 


x[n\ 



-► 


(b) 

Figure 3.13. (a) Block-diagram representation of a signal processing system; (b) 

Block-diagram representation of a discrete system. 

In general, a discrete-domain system is represented by 

y[m] = TD{x[n]} or TD{x[n]} ^ y[m], (3.23) 

where Td is a suitable mathematical operator and x[n] and y[m] represent the input 
and output sequences for such a system. Note that the sequences in general may be 
at different rates. 

2. Classification of Discrete Systems 

Discrete systems can be classihed by certain restrictions or characteristics 
placed upon, or observed concerning them. A number of such classihcations are impor¬ 
tant in signal processing, including: linear/non-linear, time-invariant/time-variant, 
causal/non-causal, stable/unstable, invertible/non-invertible systems, and systems 
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with memory/without memory. In this section, we deal with the hrst three of these 
classihcations. Others are not so important to our discussion and conform to dehni- 
tions used in many text books [Ref. 60, 61]. 

a. Linearity 

A discrete linear system is any system that satishes the generalized 
superposition principle. If the input to such a system consists of a weighted sum of 
sequences, then its output consists of the weighted sum of the system responses to 
each individual input sequence (superposition of responses). It is sufficient to dehne 
linearity in terms of the system’s response to a weighted sum of two sequences. 

Definition 20. Let xi[n] and X 2 [n] be two sequences at the same sampling rate. A 
discrete system is linear if and only if 

7D{oiiXi[n] + a2X2[n]} = ai7D{xi[n]} + a27D{x2[n]}, (3.24) 

for arbitrary input sequences Xi[n]andx 2 [n] and constants aianda 2 . 

Note that, while the output of a linear system may be at a different 
rate than the input, we make no attempt to dehne linearity in terms of two sequences 
at different rates. Indeed, the sum of two sequences at different rates has no obvious 
or unique intuitive meaning and is not necessary for our theory of multirate systems. 

Any discrete system that does not satisfy Dehnition 20 is a nonlinear 
discrete system (e.g., square law system dehned by y[n] = {x[n]y). Time dependence 
is not an issue. Both downsamplers and upsamplers (see Section III.B.1(b)) are linear 
systems. 

h. Shift-invariance 

Systems can also be characterized by the variation in their input-output 
characteristics as time evolves, and can be subdivided into shift-invariant and shift- 
dependent systems (frequently called time-invariant and time-dependent). 
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Definition 21. Let Td be a discrete system such that the input and output sequences 
are at the same rate, i.e., 7D{x[n]} {l/N}- Then the system is shift-invariant if 

and only if 

7D{x[n — iV]} ^ y[n — iV], for all integers N. (3.25) 

If the system satisfies (3.25) only for a particular value N and integer 
multiples thereof, the system is said to be periodic with period N. For systems with 
different input and output rates, shift-invariance is not defined. Periodicity, however, 
can be generalized to include these systems, as shown below. 

c. Periodic Shift-invariance 

For discrete systems with input and output signals at different rates, 
i.e., TD{x[n]} ^ y[^^], we define a particular type of shift-invariance called peri¬ 
odic shift-invariance (PSI). A system that observes this property is called a discrete 
periodically shift-invariant system. The definition of this property follows. 

Definition 22. Let 7d be a discrete system such that \i7D{x[mf\} {y[my\}, then 

the system is periodically shift-invariant if there exists integers and My such that 

7D{x[mx-Mf\} ^ y[my-My], (3.26) 

Note that when the input and output are at the same rate, a periodic 
system satisfies this definition with = My and a shift-invariant system satisfies 
this definition for all M^ such that M^ = My. The need for the more general definition 
given in (3.26) will become clear further in our development. 

d. Causality 

In defining causality, it is important to know the rates or sampling 
intervals of the sequences involved. In fact, as pointed out in [Ref. 62], “Causality is 
intrinsically related to the ordering in time [or domain] of input and output signals 
of [a] system.” Ordering is clear in a single-rate system as identical sample indices 
in each associated sequence correspond to identical points in a common domain. In 
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multirate systems, sequence sample indices must be referred back to their absolute 
scales (i.e., the discrete-domain signal) to discuss ordering. 

Definition 23. A discrete system Td is causal if and only if 


y[my\ depends only on x[mx — k] for fc = {0,1,... 


(3.27) 


where rrir = 


KyWiy 

Kx 


, and Kx,Ky are the associated decimation factors. 


A system is noncausal if it does not satisfy this definition. 

It can be seen that a discrete system is causal if and only if the discrete- 
domain signal yTy (t) at f = rOyTy depends only on values of the discrete-domain input 
signal XTx{t) for t = m^Tx < 'nXyTy. This set of input values is known as the region 
of support of the system. The concept of causality is illustrated in Figure 3.14. The 
discrete multirate system depicted has a discrete-domain input signal represented by 
the sequence x[mx] and an output signal represented by the sequence y[my]. For a 


given index in the output sequence my ^, Dehnition 23 requires that < 
for causality. 


KyCUyg 

Kx 



4m J 


W 


y[my] 


Figure 3.14. Concept of causality in a discrete multirate system comprised of a 
discrete-domain input signal x[mx] and output signal y[iTiy]. 
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3. Representation of Discrete Linear Systems 

Given a discrete linear system with input x[m] and output y[n], with, possibly, 
different sampling rates, the input-output relation or system response can be written 
as 

OO 

y[n\ = E g[n,m]x[m]. (3.28) 

m=—oo 

The term g[n,m], called the kernel or Green’s function, is the response of the system 
at point n in the output sequence to a unit impulse^ applied at point m in the input 
sequence. This formulation is a discrete-domain version of the continuous model in 
[Ref. 63, 64, 62], and the corresponding kernel has been referred to as the Green’s 
function weighting pattern response [Ref. 65]. 

a. Single-rate Systems 

When the input and output rates are the same, (3.28) can be written 
as 

CXD 

y[n] = E h[n-,m]x[n — m], (3.29) 

m=—OO 

where 

h[n-,m] = g[n,n — m], (3.30) 


is called the shift-dependent impulse response. The system is causal (see Section III.D.2(d)) 
if 


h[n] m] = 0 for m < 0. 


(3.31) 


The system is periodic (see Dehnition 21) if there exists N such that 


h[n; m] = h[n -|- iV; m] for all n. 


^The unit impulse, also know as the Kronecker delta function, is defined as 


(i[n] 


1, n = 0; 
0, n ^ 0. 


(3.32) 
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If in addition, the system is shift-invariant (see Dehnition 21), then the 
impulse response is necessarily independent of the hrst argument (n) and the output 
can be written as the familiar convolution summation 

OO CXD 

y[n]= h[m]x[n — m]= h[n — m]x[m]. (3.33) 

m=—oo m=—QC 


The system represented in (3.33) is commonly referred to as a filter because of its 
interpretation in the Fourier domain, although the term “hlter” is often used to 
apply to any system; linear or nonlinear, shift-variant or shift-invariant, in the signal 
processing literature. A hlter is said to be a finite impulse response (FIR) hlter (or 
system) if the sequence h[n] has hnite support and an infinite impulse response (HR) 
hlter (or system), otherwise. The region of support or the support of a given sequence 
is the set of values over which the sequence is non-zero [Ref. 66]. 
b. Multirate Systems 

In developing results for multirate systems, it is convenient to hrst 
describe the multirate system at the “system level” discussed in Section III.C.5 of 
this chapter. When represented at the system level, we can apply known results from 
the analysis of single-rate systems to the multirate system and then express those 
results at the signal level pertaining to the signals of interest. 

(1) System-level Representation Recall the system response 


equation (3.28) 


y[^y\ = y^j[my,mx\x[mx\, 

TUx 


(3.34) 


where y and x are at diherent rates represented by my and ma,, respectively. Also, 
recall (3.22) that a discrete-domain signal xxfit) can be represented at the system 
level 


x[mx] = x[Kxm . 


Now dehne g, the system kernel or system Green’s funetion, such that 


g[my,mx\ = g[Kymy, 



where and Ky are the decimation factors for x and y, respectively. At other values 
of it arguments, g is dehned to be 0. Now, (3.34) can be written as 

y[Kymy] = y^^glKyUiy, K^mx\ x[K^mx]. (3.35) 

rux 

If g satishes Dehnition 22 for periodic shift-invariance, then a necessary 
and sufficient condition is that 


g^Kyirriy Af^), Kxij^x A A/^,)]. 

where and My are the number of samples per system period for x and y, respec¬ 
tively. Then 


g[Ky{my -h My), Kx{mx + il4)] = g[Kymy -F N, K^rrix + A^] (3.36) 

where we recall (3.17) that the system sample period is given by = KyMy = K^M^', 
therefore, 'g has period N in both arguments. 

Now, we dehne 

h[Kymy-, K^rrix] = g[Kymy, Kyruy - K^rrix], (3.37) 

or, equivalently, 

g[Kymy-, K^rrix] = h[Kymy, KyfUy - K^rrix], 

where h is called the system impulse response. Notice that since g satishes (3.36), 
then 

h[n-, 1] = 'g[n, n — 1] = 'g[n + N,n — I + N] = h[n -|- A^; /], 

thus /i[n; /] is periodic in its hrst argument only. Note, also, that if the system is 
causal, then h[n; 1] must be 0 for / < 0. Now we can write 

g[Kymy-, K^mx] = h[Kymy + N] Kymy - Kxmx], (3.38) 

where N is the system sample period. 
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Now, (3.35) can be written as 

ylKyUiy] = ^ h[Kymy] Kyiriy - (3.39) 

rux 

Since h is periodic in iV, we can define 

= /i[n; /] where p = (n) (3.40) 


and N is the system sample period. Substituting into (3.39) yields the general filter 
equation 

IX] “ K^tUx] ^Kxrrix ], Ky\n, 


y[n] = 


1 


(3.41) 


undefined, otherwise 

where p = (n) 

If we desire a causal filter (see Section III.D.2(d)), then 

KyUly 

rrixTx < rUyTy or nix < —^—, 


(3.42) 


and 


( I I 

L if,,. J 

[n - KxUix] x[Kxmx ], 

mx=—oo 



Ky\n, 


(3.43) 


y undefined, otherwise 

where p = (n) The upper limit on the sum is insured if h [1] = 0 for I < 0. Notice 
that when Tx = Ty ^ Kx = Ky (single-rate system), (3.41) and (3.43) simplify to the 
expected convolution sums. 


Example 11. Consider Figure 3.15, where two signals are represented on their as¬ 
sociated grids. The upper grid (a) represents a signal y that has an associated dec¬ 
imation factor Ky = 3, and the lower grid (b) represents a signal x that has an 
associated decimation factor Kx = 2. Recall (3.18) that the system period is given by 
N = \cm{Kx, Ky); therefore, N = 6. 
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nT„ 
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fiT 


0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 

(c) 


Figure 3.15. (a) Discrete-time signal y[n] with decimation factor Ky = 3; (h) 

Discrete-time signal x[n] with decimation factor =2; (c) System grid. 


If we desire the estimate for y at my = A, using a causal filter, where 
the filter order P = 3, then we can write (3.43) as 


LfJ=6_ 

^[12] = h^^\l2 — 2mx]x[2mx], where p = (12)g = 0. 

mx=3 

therefore 

y[12] = 7):^°^[4]a;[8] +h^''\2]^10] + 7):^°^[0] a;[12], 


Notice that the linear combination is in terms of the system-rate parameters (n is the 
system time index). Once y[12] is computed, recall that y[my] = y[myKy]; therefore, 
y[A]=y[12]. 
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If we desire the estimate for y at my = 5, 

LfJ=7_ 

y[l5] = ^ — 2mx]x[2mx], where p = (15)g = 3. 

mx=S 

therefore 

y [15] = h^^^ [5]x[10] + h^^^ [3]x[12] + h^^^ [l]x[14]. ■ 

(2) Signal-level Representation Often, it is desirable or more 
convenient to deal in terms of the actual signal parameters rather than system pa¬ 
rameters; therefore, we develop (3.41) and (3.43) in terms of the individual signal 
parameters. To do so, we introduce the following substitutions and change of vari¬ 
ables. 

First let us examine the system impulse response h, 
h^^^ [Kymy - Kxmx] = h[Kymy -F N; Kymy - Kxmx]- 
Consider the first argument, where Ky has been factored out and we recall (3.17) 

Kymy + N = KyirUy + N/Ky) = KyirUy + My). 

Since this argument is periodic. 


Kyirriy -\- N/Ky) <—>■ Kyiniy) ^ 


Now, consider the second argument, where Kx has been factored out 

' K.m,^ 


Kyiriy KxTtIx Kx 


^yliuy 


Kx 


— m. 


Recalling that 


and 


y[my] = y[myKy], x[mx\ = x[mxKx], 
h[Kymy] Kxmx\ = h[my-,mx], 
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we can write the impulse response as 




KyUly 

K. 



and we can write 


where l = {my)^ 


y[my] = 


KyUly 


— m. 


x[mx]i where I = (my) 




If we introduce the following change of variables 


k = 


KyUly 

Kx 


rrin 


(3.44) 


we can write the general signal-rate hltering equation as 


= X] h^^\mx]x 

rrix 

Kymy 
[ Kx \ 

- nix 

, where/ = 


and the signal-rate causal filtering equation as 


1 _ 

H 

Kymy 

Kx 

-rux 

, where/ = 

mx=Q L 


- 



(3.45) 


(3.46) 


The lower limit {mx = 0) on the summation is equivalent to the causal condition 
h^^^rUx] = 0 for rrix < 0. Notice that My, the number of signal samples (y) per 
system period, is the period of the hlters required to form y[my], referred to as the 
cy do stationary period in [Ref. 9]. Additionally, note that when Tx = Ty, then 
Kx = Ky and (3.45) and (3.46) simplify to the expected convolution sums for a 
single-rate system. 


Example 12. Let us illustrate with the same example previously discussed and illus¬ 
trated in Figure 3.15. Recall that signals y and x have associated decimation factors 
Ky = 3 and Kx = 2, respectively, and that the system period N = 6. Further, recall 

, .X N 6 

that My = = - = 2 . 

JV 71 3 
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If we desire the estimate for y at my = 4, using a causal FIR filter, 
where the filter order P = 3, then we can write (3.45) as 

2 

y[^ = h^^\mx\x[I) — mx], where I = = 0, 

mx=0 

therefore 

y[4] = /i(o)[0]a;[6] + h^^\l]x[5] + h^^\2]x[4\. 

If we desire the estimate for y at my = 5, 

2 

l/[5] = h^^\mx]x\l — mx\, where I = ( 5)2 = 1 , 

mx=0 

therefore 

y[5] = /i^^^[0]a;[7] + [l]a;[6] + h^^'^[2]x[5]. ■ 

E. MATRIX REPRESENTATION 

In the analysis of multirate systems, linear algebra concepts provide a useful 
framework to represent or model basic multirate operations, such as downsampling, 
upsampling and linear hltering. This section develops a systematic methodology, 
which is used in other sections to further represent multirate systems, an extension 
of work presented by [Ref. 67]. 

1. Decimation 

Decimation is the process of digitally reducing the sampling rate of a signal 
[Ref. 65], and the basic element used in digital systems is the downsampler or dec- 
imator, shown in Figure 3.16. This element extracts every sample of its input 
and, as a result, decreases the sampling rate by Its operation can be expressed 
mathematically as 

y[n] = x[Mn] (3.47) 
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x[n] 


Figure 3.16. M-fold downsampler. 


where n is an integer. 

If the input and output signals are represented as vectors, then this operation 
can be expressed as a linear transformation of the form 

y = 0-{M,iV,x} (3.48) 

where T{-} represents an arbitrary linear operator that is dependent on the decimation 
factor M and the order of its associated vectors. If the downsampler takes input vector 
X G where 

X = [ a;[0] a:[l] ... x[M] x[M + 1] ... a;[iVM — 1] Y 

and produces an output vector y G where 

y=[a;[0] x[M] x[2M] ... x[(iV - 1)M] ]^, 

then the operator can be expressed as 

7 =-Dm, (3.49) 

and the transformation can be written as 

y = Dmx. (3.50) 

The constant matrix Dm G is called a decimation matrix and is dehned in 

terms of a Kronecker product (introduced in Section II.B.2) as 

Dm = 1® (3.51) 

where I is the N x N identity matrix and is an M x 1 index vector with a 1 in the 
hrst position and O’s elsewhere. 
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As an example, consider the case where M = 3 and x G The decimation 
matrix D 3 G is given by 


Ds 


1 0 0 
0 0 0 
0 0 0 
0 0 0 


0 0 0 
1 0 0 
0 0 0 
0 0 0 


0 0 0 
0 0 0 
1 0 0 
0 0 0 


0 0 0 
0 0 0 
0 0 0 
1 0 0 


which resnlts in y G 


2. Expansion 

Expansion and decimation are dnal operations. Expansion is the process of 
digitally increasing the sampling rate of a signal, and the basic element nsed is the 
npsampler or expander, shown in Fig 3.17. This element inserts L — 1 zeros after every 
sample of the inpnt and, as a resnlt, increases the sampling rate by L. Mathematically 
the process is described by 


y[n] 



if n is an integer mnltiple of L; 
otherwise. 


(3.52) 


or 

CXD 

y[n] = E x{k)6{n — kL). (3.53) 

k=—oo 


x[n] 



y[n] 


► 


Fignre 3.17. L-fold expander. 

If the inpnt and ontpnt signals are represented as vectors, then this operation 
can be expressed as a linear transformation of the form 


y = T{L,iV,x} 


(3.54) 
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where T{-} represents an arbitrary linear operator that is dependent on the upsam¬ 
pling factor L and the order of its associated vectors. If the expander takes input 
vector X G where 


= a;[0] a;[l] ... x[N — 1] ]' 


and produces an output vector y G where 

L-l L-l 

y = I a;[0] 0 0 a;[l] 0 0 x[2] 


L-l 

x[{N - 1)] o'... 0 


then the operator can be expressed as 


T=Ul, (3.55) 

and the transformation can be written as 

y = Ulx. (3.56) 

The constant matrix Ul G is called an expansion matrix and is dehned in 

terms of a Kronecker product as 


UL = I®t (3.57) 

where I is the N x N identity matrix and t is an L x 1 index vector with a 1 in the 
hrst position and O’s elsewhere. 

As an example, consider the case where L = 3 and x G The expansion 
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matrix U 3 G is given by 


10 0 0 
0 0 0 0 
0 0 0 0 
0 10 0 
0 0 0 0 
0 0 0 0 

U3 = 

0 0 10 
0 0 0 0 
0 0 0 0 
0 0 0 1 
0 0 0 0 
0 0 0 0 

which results in y G 

Notice that, if L = M, the expansion matrix is related to the decimation 
matrix through matrix transposition as 

Dm = {Uif for M = L- (3.58) 

in other words, the operation of decimation and expansion are complementary pro¬ 
cesses when M = L. 

3. Sample Rate Conversion with Delay 

In the analysis of various multirate implementations, notably those involving 
polyphase decompositions (efficient implementations) [Ref. 68], signal delay must be 
incorporated in the system model and forms another basic building block of multirate 
systems as shown in Figure 3.18. In discrete systems, a signal delay^ is simply a shift 

^Despite the name, there is no restriction in considering advancement in sample index. This will 
be clear by the sign associated with the delay. 
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of the associated sequence to the right or left by some integer number of samples M 
and is represented by the usual 

x[n — M] <->■ z~^ 

and 

x[n + M] <—> 

where z~^ indicates a delay or shift to the right, and z^^ indicates a negative delay 
or shift to the left. 

In conjunction with signal delay, the sample rate conversion transformations 
of (3.50) and (3.56) can be generalized. Given a delay A; G {0,1,..., M — 1}, then 

y = Dm X. (3.59) 

The constant matrix G is called an decimation matrix with delay and is 

defined in terms of a Kronecker product as 

= I 0 4 (3.60) 


where I is the N x N identity matrix and is an M x 1 index vector with a 1 in the 
{k + 1)*^ position and O’s elsewhere. 


x[«] 



♦ 



y[«] 


► 


Figure 3.18. M-fold decimator with delay. 


In a similar manner, 

y = (3.61) 

The constant matrix G is called an expansion matrix with delay and is 

defined in terms of a Kronecker product as 

=l®Lk (3.62) 
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where I is the N x N identity matrix and is an L x 1 index vector with a 1 in the 
{k + 1)*^ position and O’s elsewhere. 

Consideration of (3.60) and (3.62) shows the duality of the two operators, for 
M = L, and results in the general expression 

= = l®Lk forA;e{0,l,...,M-l}. (3.63) 

4. Linear Filtering 

It is also useful to provide a matrix representation of linear hltering. Recall 
from (3.33) that a P length causal FIR hlter is fully specihed by its P hlter coefficients 
(impulse response). If we denote these coefficients h G as 

hp = [ h[0] h[l] ... h[P-l] f 

and the related input data sequence xp[n] G as 

xp[n] = [ x[n — P + 1] ... x[n — 1] x[n\ , 

then the output of the hlter at index n, y[n\ G M can be represented by 

y[n] = hpXp[n] = hpXp[n]. (3.64) 


F. CHAPTER SUMMARY 

This chapter establishes the theory of multirate systems and provides the foun¬ 
dation upon which the remainder of this work is built. The concept of a multirate 
system is introduced, with various practical examples. A multirate system is formally 
dehned, the notion of rate is discussed, and the basic multirate operations of down- 
sampling and upsampling are introduced. A conceptual framework for the analysis 
of multirate systems is developed, which enables a systematic extension of optimal 
estimation and linear hltering theory to multirate systems. The characterizations of 
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constituent multirate signals introduced include system rate, decimation factor, sys¬ 
tem period, and maximally-decimated signal sets. Further, this chapter introduces 
the concept of a system grid, which allows representation of the various signals in a 
multirate system on a common domain. 

The system theory of multirate systems is also developed, including the con¬ 
cepts of linearity, shift-invariance, period shift-invariance, and causality. In addition, 
the input-output relation of a multirate system is discussed in terms of the system re¬ 
sponse and its associated Green’s function and is adapted to the multirate problem in 
both system-level and signal-level representations. Further, the relationship between 
linear periodically time-varying hlters and their linear time-invariant equivalents is 
discussed. 

Finally, the basic multirate operations are analyzed from a linear algebraic 
point of view, and matrix representations for the operations of downsampling, up¬ 
sampling, sample rate conversion and linear hltering are presented. 
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IV. MULTIRATE OPTIMAL ESTIMATION 


Signal or image reconstruction can be viewed as a problem in signal estima¬ 
tion, where a related low-rate (low-resolution) signal or signals is used to estimate 
an underlying high-rate (high-resolution) signal. From this perspective, the observa¬ 
tion signal or signals, and desired signal form a multirate system and the theory of 
multirate systems developed in Chapter III can be used to extend single-rate signal 
estimation theory to the multirate case, which is the concern of this chapter. 

A. SIGNAL ESTIMATION 

“Estimation is the process of inferring the value of a quantity of interest from 
[typically] indirect, inaccurate, and uncertain observations [Ref. 69].” A pictorial 
representation of this concept is depicted in Figure 4.1, where the source d represents 
some quantity of interest (unknown parameter, random variable, random process, 
etc.), X are the observations, related to d, and (i(x) is the estimate. Note that, since 
X and sometimes d are considered to be random variables, d, which is a function 
of the observations, is also a random variable. In the field of signal processing, 
these quantities of interest are signals, and a major emphasis of research is on signal 
estimation, where one signal is estimated from some other related signal or signals. 
The desired signal may be corrupted by distortion or interference and is usually 
unobservable (at least at the moment when the estimate is desired). In [Ref. 70], a 
number of typical signal estimation applications are presented including: the recovery 
of a transmitted signal from a distorted received signal, subject to amplitude and 
phase distortions and additive white noise over the communications channel; and 
image restoration of an image recorded by an imaging system that introduces blurring, 
nonlinear geometric distortions, and additive white noise. 
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Figure 4.1. Concept of estimation. 


B. OPTIMAL FILTERING 

When the quantity of interest is a random signal, optimal hltering provides 
a framework for signal estimation. Optimal filtering is an area of signal processing 
study that is concerned with the design of filters to process a class of signals with 
statistically similar characteristics [Ref. 5] that are “best” in some sense, in terms of 
stated optimality criteria (i.e., minimum mean-squared error). The optimal filtering 
problem is posed in the following manner. Suppose that a random process x[n\ is 
observed, which is related to another random process d[n\ that cannot be observed 
directly. A general expression for the estimate is 

d[n] = 0[{a:[n]}]. (4.1) 

The optimal filtering problem is concerned with finding the appropriate functional 
0[-] that provides the best estimate of the desired signal d[n]. When 0[-] is linear. 
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the functional is commonly referred to as a (linear) filter. This concept is depicted in 
Figure 4.2. 



Figure 4.2. General single-rate optimal filtering problem. When (()[■] is linear, the 
functional is commonly referred to as a linear filter. 


Typically, these hlters are constrained to be a linear function of the observa¬ 
tions and hence can be written in the form (recall (3.28)) 

OO 

d[n] = g[n]m]x[m], (4.2) 

m=—oo 

where the sequence g[n]m], which may be finite or infinite, is the familiar Green’s 
function. This class of optimal hltering is called optimal linear filtering. 

Often, the optimality criterion is based on minimization of the mean-square 
error (MSE) between the desired and estimated signals, i.e., 

a2[u] = £{|£[n]|2}, (4.3) 

where 

e[n] = d[n] — d[n] (4.4) 

is the error in estimation at the sample. When coupled with this MSE optimal¬ 
ity criterion, optimal linear hltering is commonly referred to as Wiener filtering in 
recognition of the pioneering work of Norbert Wiener [Ref. 71, 72] in the 1940’s.^ Of 
great signihcance is that only second-order statistics are required in determination 

^Wiener’s pioneering work addressed the optimal filtering problem in continuous time, but is 
easily translated to discrete time, as is common now. Kolmogorov, in Russia, worked on the discrete¬ 
time problem for time series [Ref. 73] and appears to have preceded (or at least matched) Wiener 
in formulating and solving the problem (for discrete time). 
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of these optimal filters [Ref. 61]. Further, under the condition of joint wide-sense 
stationarity (JWSS) between the input x[n] and the desired process d[n], the hlter 
is shift-invariant (see Section III.D.2(b)). For our purposes, we develop the general 
discrete FIR Wiener hltering equations, then the hltering equation for JWSS observa¬ 
tion and estimate signals. In order to facilitate such developments, the orthogonality 
principle [Ref. 5] is stated. 

1. Orthogonality Principle 

Theorem 3 (Orthogonality Principle). Let e[n] = x[n]— x[n] be the estimation 
error. Then, 

1. the optimal linear hlter with coefficients h[0], h[l],..., h[P — 1] minimizes the 
error variance if the hlter coefficients are chosen such that E{e[n]x*[n — i]} = 
8,{x[n — f]£*[n]} = 0, i = 0,1,..., P — 1, that is, if the error is orthogonal to 
the observations. 

2. The minimum error variance is given by = £{e[n](J*[n]} = £{(J[n]e*[n]}. 
The proof of this result can be found in several places [Ref. 5, 60, 74, 70]. 

2. Discrete Wiener Filter Equations 

The equations whose solution provides the optimal hlter are known as the 
Wiener-Hopf equations. In developing the “ordinary” (single-channel/single-rate) 
Wiener-Hopf equations, let us dehne the estimate in terms of a linear FIR hlter that 
operates on x[n] and the P — 1 previous values of the process. The estimate (from 
(3.29)) is given by 

p-i 

d[n] =h[n;m]x[n — m], (4.5) 

m=0 

where h[n]m] is the shift-dependent impulse response sequence of the optimal hlter 
and P is the associated hlter order (length of the impulse response sequence). If 
x[n — i] for i = 0,1,. .., P — 1 represents any one of the observations of the sequence 
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x[n], then by applying the orthogonality principle (Theorem 3), we can write 


p-i 


E{e[n]x*[n — i]} = E, < i d[n] —h[n] m]x[n — m] x*[n — i\\ = Q 


m=0 


or 


p-i 


R^[n — m',i — m] — i 


i = 0,1. 


(4.6) 


(4.7) 


m=0 


This equation is the desired discrete Wiener-Hopf equation and can be written in 
matrix form. For example, for P = 3, we have 

Rx[n-,0] Rx[n-l-,-l] Rx[n-2-,-2] h[n;0] Rdx[n-,0] 

Rx[n-,l] Rx[n-1;0] Rx[n-2]-l] h[n; 1] = Rdx[n]l] ■ (4.8) 

Rx [n; 2] [n - 1; 1] Rx [n - 2; 0] /i[n; 2] Rdx [n; 2] 

By applying the second part of the orthogonality principle, the expression for 
the minimum mean-square error can be found as 


p-i 


(T^min W ~ E{e[n]d*[n]} = £ < j d[n] — /i[n; m]x[n — m] d*[n 


(4.9) 


m=0 


or 

p-i 

^Ln W = Mn] 0] - h[n] m]Rax[n] m]. (4.10) 

m=0 

If the random processes are JWSS, then the statistical properties of the as¬ 
sociated system do not change with n, the estimate can be generated by a linear 
shift-invariant FIR hlter, and the minimum mean-square error becomes constant. 

The matrix form of the Wiener-Hopf equation (4.7) can be derived directly if 
(4.5) is written in vector form. 


where 


d[n] = (x[n])^h. 


x[n 


x[n — P -I- 1] 
x[n — P] 


x[n 


(4.11) 


(4.12) 
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with x[n] as defined in Section II.B.3, and 

fc|0] 

Ml] 

h\p - 1] 

or in the general non-stationary case 

/i[n; 0] 

/i[n; 1] 

h[n; P — 1] 

Again, evoking the orthogonality principle, we can write 

R*h = Frfx, (4.15) 

where 

R:e = £{x[n]x*^[n]}, (4.16) 

and 

rdx = £{d[n](x[n])*}. (4.17) 

Equation (4.15) has the specific form (4.8) (for P = 3). The minimum mean-square 

error can be written as 

^Ln W = [n; 0] - f . (4.18) 

C. MULTIRATE OPTIMAL FILTERING 

1. Single-channel, Multirate Estimation Problem 

In the single-channel, multirate optimal filtering problem, the observation sig¬ 
nal and the desired signal are at different rates and our goal is to form an estimate 


(4.14) 



(4.13) 
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) x[m] 

Filter 



(^[■] 



d[n\ 


Figure 4.3. General single-channel, multirate optimal filtering problem. Note that 
the estimate and observation signals may be at different rates. 


of the underlying desired signal from the observation signal. This notion is depicted 
in Figure 4.3 where x[m] represents the observation signal at rate d[n\ represents 
the desired signal at some other rate Fd, and d[n\ represents the estimate. A general 
expression for the estimate is 

d[n] = 0[{x[m]}], (4.19) 

where the scalar estimate is a function of the sequence x[m]. Again, we seek the 
appropriate functional that provides the best estimate of d[n\. When 0[-] is linear, 
the functional is, again, referred to as a linear hlter. 

In this work, we have developed two approaches to formulating the required 
estimate. The hrst approach is based on a causal FIR Wiener hltering model, the 
second on a non-causal FIR Wiener hltering model. The latter is of particular interest 
since it is the basis for the two-dimensional work that follows (see Chapter V). 

a. Index Mapping 

The FIR linear hltering problem involves computing the linear combi¬ 
nation of some hnite number of observations to determine an estimate of a desired 
signal at a particular sample index, n = uq. 

In the single-rate case, the region of support of the hlter (see Sec¬ 
tion III.D.3) is unambiguous. In the case of a causal hlter, the hlter uses x[n(^ 
through x[nQ — P -|- 1] to estimate (i[no]. In the multirate case, determining this re¬ 
gion of support requires a way of relating corresponding sample indices between the 
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desired and observation seqnences. Developing this methodology, referred to as index 
mapping, is the concern of this section. 

Let ns motivate the discnssion by considering Fignres 4.4 and 4.5, where 


d\n\ 
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x\n\ 
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Fignre 4.4. An illnstration of ordinary cansal FIR Wiener hltering and the relationship 
between samples of seqnences d[n] and x[n], P = ?>. 


the horizontal axes represent time. The hrst hgnre depicts single-rate estimation of the 



x[m\ 
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Fignre 4.5. An illnstration of single-channel, mnltirate cansal FIR Wiener filtering 
and the relationship between samples of seqnences d[n] and x[m], P = 2. 


desired seqnence d[n] at sample index n = k with a cansal FIR filter of order P = 3. 
Since both seqnences are at the same rate, the region of snpport, the region encom¬ 
passed by the rectangnlar “filter”, corresponds to observations x[k], x[k — 1], x[k — 2]. 
If we consider the nnderlying discrete-domain signals, we notice that both signals are 
defined on the same “signal domain”, Since both signals are dehned on 
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the same domain, the index value k for either sequence represents the same point in 
time. 

Figure 4.5 depicts multirate estimation of the desired sequence d[n] at 
sample index n = k with a causal FIR hlter of order P = 2. Here, we notice that the 
underlying discrete-domain signals are not dehned on the same domain. So, consider 
a particular index value k. In general, the time at which d[k] occurs {kTd) is not the 
same as the time at which x[k] occurs {kT^). There may be some other index value 
I (see Figure 4.5), however, such that d[k] and x[l] occur at the same time or very 
nearly at the same time. In an optimal hltering problem, we would generally like 
d[k] and x[l] to be as close as possible. In addition, we may also require the hlter to 
be causal (see Dehnition 23). Establishing this correspondence is referred to as the 
index mapping problem. 

Let Trf = {nTd] n G Z} represent the signal domain that corresponds to 
an arbitrary discrete-domain signal with sampling interval Td, and Ta, = {mTx]m G 
Z} represent the signal domain that corresponds to its associated observation signal 
with sampling interval T^. Also, dehne a distance metric [Ref. 54] D[n,m] called the 
index metric as 

D[n,m] = \nTd-mTa:\ = T\Kdn - K^m\, (4.20) 

where T is the system sampling interval. The index metric gives us a notion of 
distance between any two indices from different signal domains. In fact, \Kdn — K^ml 
is the number of system samples between such indices. Figure 4.6 illustrates this 
concept of distance for two arbitrary sample indices, no and mo. A normalized plot 
of the index metric, for = 3 and Kd = 1, is depicted in Figure 4.7(a). In this 
case, D[n,m] = \n — 3m|, for 0 < n,m < 30. Each line contains the locus of points 
associated with a particular value of n. Also, note that the index metric is identically 
zero D[n, m] = 0 when n = 3m. 

The index mapping problem is now stated. Given a discrete-domain 
signal d at rate Fd and its associated observation signal a; at a different rate Fx, 


81 



Td 






( ( 






0 4 

d 2 

Td 


) ) 

• • • 



V 

^d 




( ( 

1 

1 



) ; 

1 


0 1 

• • • 

X 

! “oTk 


I 

I 

I' 


D{n^,m^) = T\K^n^-K^m^ 
Figure 4.6. Notion of distance between indices no and mo. 


for a given sample index n = no, associated with d, find the observation sample 
index m = mo, associated with x such that the related index distance D[no,mo] is 
minimized, subject to certain constraints. Intuitively, this makes sense since the best 
estimate of d[no] should typically result from the closest observation samples. 

If the problem of interest involves causal filtering (see Definition 23), 
then the required constraint on the minimization is 

K 

nTd > mTx or n > (4-21) 

Kd 

This states that (on the system level) the output must not precede the observations. 
Therefore, the problem is posed as 


mmD{n, m) 


mKx 

subject to n > ———. 

Kd 


(4.22) 


The solution to this minimization is expressed as a mapping from the 
estimation index to the observation index as 


^^causal • ^ 


nKd 


or 


^ -^causal (^), 


(4.23) 


(4.24) 
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where Kd and are the associated decimation factors. 

This minimization is illustrated in Figure 4.7(b), which displays a plot 
of the index metric, again for = 3 and Kd = 1, where the value of n has been 
specihed as n = 5. In this case, the causality constraint requires that n > 3m, which 
is indicated by the shaded region in the lower hgure. Thus, the minimization pertains 
only to that portion of D[5,m] contained within this region. By inspection, we see 
that when m = 1, the index metric D[5,m] = 2 is at its minimum. 

Continuing this example, for = 3 and Kd = 1, with the causality 
constraint, we have the following mappings from the set of estimate signal indices to 
the observation signal index as shown in Table 4.1. 



Table 4.1. Causal mapping from a set of estimate signal indices to the associated 
observation signal index. 


If there is no causality constraint (as in smoothing or image processing), 
then the minimization is unconstrained and the problem is posed as 


min D[n, m . 


(4.25) 


The solution to this minimization is expressed as a mapping from the 
estimate signal index to the respective observation signal index as 


: n 


{n+m)Kd 

K^ 


or 


m = Mnc(n), 

where Kd and K^ are the associated decimation factors. 


(4.26) 


(4.27) 
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Figure 4.7. (a) Normalized plot of D[n, m] in 3 dimensions, (b) Plot of D[n, m] versus 
m for n = 5. 
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Again, from Figure 4.7(b), we can illustrate this minimization, again at 
n = 5, but now with no causality constraint. In this case, when the observation signal 
index is m = 2, the index metric is D[5, 1] = 1, and the index metric is a minimum. 
Note that this value of the observation signal index m is different from the previous 
solution (a consequence of the causality constraint in the former problem). 

Continuing this example, for = 3 and Kd = 1, without the causality 
constraint, we have the following mappings from the set of estimate signal indices to 
the observation signal index as shown in Table 4.2. 



Table 4.2. Non-causal mapping from a set of estimate signal indices to the associated 
observation signal index. 


b. Single-channel, Multirate Wiener-Hopf Equations 

Recall that the input-output relationship for the general multirate £1- 
tering problem is given by (3.45) and is repeated here for convenience 


y[my\ = ^ ^ 




— m. 


where I = (my) 




This equation can be written in terms of the single-channel, multirate estimation 
problem as 


d[n] = ^ h^''\m]x 

m 


- KdU 

1 K, . 



where I = (n) 


We can further generalize the form of the linear multirate estimate if we consider the 
mapping issues discussed in Section IV.C. 1(a). We can write (3.45) in terms of an 
arbitrary mapping function as 


d[n] = E 


m 


(4.28) 
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where M[-] is the appropriate mapping function (causal, non-causal), and I = 

With a hnite variant of (4.28), in the usual manner, we evoke the 
orthogonality principle. Thus, for observations a;[M[n] — i] for f = 0,1,..., P — 1, we 
can write 


8,{e[n]x* M[n] —i} 

= £ / (d[n\ — ^ h^''\m]x M[n] — m 


p-i 


m=0 



for i = 0,1,... ,P-1, (4.29) 


or 

p-i 

^^Rx[y^[n]—m]i — m]h^’'\m] = Rdx[n]n — 'M.[n]+i]] f = 0,1,..., P — 1, (4.30) 

m=0 

where I = (n) This is called the discrete single-channel, multirate Wiener-Hopf 
equation. 

Applying the second part of the orthogonality principle, we hnd the 
following expression for the minimum mean-square error 

p-i 

^Ln W = ^d[n-, 0] - h(') [m]R*^Jpi- n - M[n] + m]. (4.31) 

m=0 

Notice that if the observation sequence is stationary, then the associated 
correlation function is independent of the mapping function and the Wiener-Hopf 
equation can be written 

p-i 

Rx[i — = Rdx[n', n — M[n] -|- m]; f = 0,1,..., P — 1, (4.32) 

m=0 

where I = {n) 

The Wiener-Hopf equation expressed in its matrix form, for P = 3, 
when the observation process is stationary: 


Rx[0] Rx[-l] Rx[-2] 


h«[0] 


Rdx[n]n- M[n]] 

Px[l] Pz[0] Pz[-1] 


h«[l] 

= 

Rdx[n-,n- M[n] -h 1] 

Rx[2] Rx[l] Rx[0] 


h«[2] 


Rdx[n-,n- M[n] 2] 


(4.33) 



where I = (n )Notice that the correlation matrix is Toeplitz. 

As in the ordinary Wiener hltering problem, the Wiener-Hopf equation 
can be derived in matrix form by writing (4.28) as 


d[n\ = (x M[n] 


(4.34) 


where if 


then 


X M[n] 


a;[M[n] — P + l] 
a;[M[n] — P] 


X M[n] 


X M[n] 


a;[M[n] — l] 




X M[n] — P + l] 


with x[M[n]] as dehned in Section II.B.3, and where is dehned as 


(4.35) 


(4.36) 


h(0 


h«[0] 

h«[l] 


/i(') [P - 1] 


(4.37) 


where I = {n) 

Again, evoking the orthogonality principle (Theorem 3), the Wiener- 
Hopf equation (4.30) can be expressed in its matrix form as 


where I = {n )and 




dx 3Vt[n] ’ 



£{(i[n]x* M[n] }. 


(4.38) 


(4.39) 
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Equation (4.38) has the specihc form (4.33) (for P = 3). The minimum 
mean-square error can be written as 


In] = Rd\n: Ol — r 

Emin I- J “ L 5 J 


dx 


M[r 


hh)-* 


(4.40) 


where I = {n) 

c. Matrix Approach to the Single-channel, Multirate 
Wiener-Hopf Equations 

Another approach in developing the multirate Wiener-Hopf equations 

involves the matrix representation concepts of Section III.E and results presented 

in [Ref. 9]. Again, consider a multirate system comprised of an observation signal 

X at rate and a desired signal d at rate Fd, and let the order of the associated 

hlter be denoted P. Recall from (3.10), (3.14) and (3.18) that the system rate is 

P = lcm{Fx, Fd), the associated decimation factors are Ki = F/Fi, and the system 

sample period is iV = lcm{Kx, Kd). The period of the hlters required to form the 

estimate d (see Section III.D.3.b(2), cyclostationary period) is K, and in the case of 

N 

single-channel, multirate estimation, K = 

Kd 

For any index n, the observation signal x, represented by the sequence 
x[n], can be expressed 







(4.41) 


where 


X n = 


x[n\ 

x[n — 1] 
x[n — PKx + 1] 


and is a decimation matrix with delay (3.60). In this context, the decimation 
matrix with delay is a mapping or transformation that extracts the samples in the 
required causal region of support from the observation vector. 



The optimal estimate is formed from a linear combination of the ap¬ 
propriate observation samples and associated hlter coefficients and is given by 

dk[n] = k = {0,l,...,K-l}, (4.42) 

where the reversal operation x[n] is dehned in Section II.B.3. Note that in this 
formulation, all signals and computations are at the system rate, and, as a result, 
we are solving a more general problem since we estimate [n] at every point on the 
system grid not just on the signal domain The desired estimate is recovered 
by decimating the result by Kd- 

In the usual manner, we form the error in estimation and evoke the 
orthogonality principle. The Wiener-Hopf equation (4.30) can be expressed in its 
matrix form as 

A; = {0,1,..., ir - 1}, (4.43) 

or in terms of single-rate statistical parameters (from (4.16) and (4.17)), 

A; = {0,1,..., 77 - 1}. (4.44) 

Applying the second part of the orthogonality principle yields an ex¬ 
pression for the minimum mean-square error, 

al = i?40] - f A; = {0,1,..., 77 - 1}, (4.45) 

where the subscript k on the mean-square error reflects the periodically time-varying 
nature of the error. 

2. Multi-channel, Multirate Estimation Problem 

In the multi-channel, multirate optimal hltering problem, there are multiple 
observation signals and these signals are allowed to be at rates different from the 
desired signal. The goal is to estimate this underlying signal from the set of related 
multirate observation signals. This notion is depicted in Figure 4.8 where there are 



M observation signals at rates denoted by m*, which may all be different. A general 
expression for the estimate is 

d[n] = (f) [{xi[mi];i = 0,1, • • •, Af - 1}], (4.46) 

where the scalar estimate is now a function of a set of multirate observation signals. 
When the functional 0[-] is linear, then the estimate is formed from linear filtering. 



Figure 4.8. General multirate optimal hltering problem with M multirate observation 
signals. 

a. Multi-channel Index Mapping 

In this section, we develop the mapping relationship between the esti¬ 
mate and observation sample indices in a multi-channel, multirate system in order to 
determine the required regions of support in linear filtering. In doing so, we develop 
a set of mapping functions for the causal and noncausal FIR filtering problems. In 
this context, the index mapping problem is, for a given index n, associated with the 
estimate, to find the M observation indices m*, such that the related index distances 
D[n, rrii] are minimized, subject to certain constraints, where i = 0,1,..., M — 1. 

If the problem of interest involves causal hltering, then the required 
constraint on the minimization is 

K 

nTd > rriiTi or n> — -nii, for i = 0,1,.. ., M — 1. (4.47) 

Kd 

Therefore, the problem is posed as 

Ki 

mmD{n,mi) subject to n >-^m*, for i = 0,1,..., M — 1. (4.48) 

rrii Kd 
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The solution to this minimization is expressed as a mapping from the 


estimate index to the observation indices as 


'^^causal • ^ 


nKd 

Ki 


(4.49) 


or 

rrii = Mcausai(n), (4.50) 

where Kd and Ki are the associated decimation factors, and z = 0,— 1. 

The multi-channel, multirate causal hltering problem is depicted in 

Figure 4.9. 

d\n\ 





xjwj 



m. 


Figure 4.9. Concept of index mapping in multi-channel, multirate FIR Wiener filter¬ 
ing. 


If the problem of interest involves noncausal FIR filtering, then there 
is no constraint on the minimization, and the problem is posed as 


minD[n^mi] for f = 0,1,... , M — 1. (4-51) 

mi 
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By examining the appropriate plots of -D[n, m*], the solution to this 
minimization is found as a mapping from the high-rate index to the respective low- 
rate indices as 

(4. 

or 





: n 


{n+m)K, 

Ki 


rrii = Mnc(n), (4.53) 

where Kd and Ki are the associated decimation factors, and i = 0,— 1. 

h. Multi-channel, Multirate FIR Wiener Filtering model 
Recall that, for the single-channel FIR Wiener hltering model (4.28), 
the estimate is written as 

p-i 

d[n] = ^ h^’'\mx]x [M[n] — tUx] , 

m ^=0 

where I = {n)ji_, and M[n] is the required mapping function (causal, non-causal, 
etc.). 

In the multi-channel model, we extend or generalize (4.28) to include 
multiple observations as 

M-l P-1 

d[n] = ^ ^ h^^\m]xi [Mj[n] — m] , (4.54) 

i=0 m=0 

where p = (n) ji_, P is the hlter order, and Mj[n] is the mapping function associated 

Kd 

with the channel. 


c. Multi-channel, Multirate Wiener-Hopf Equations 

Using equation (4.54), we form the error in estimation and evoke the 
orthogonality principle. Thus, we can write 


t{e[n]x* Mj[n]—j } 


M-l P-l 


£ d[n] [sjxr Mr[?^] — s I a;* M* [n] — j 


r=0 s=0 


= 0 , 


0<i<M-l,0<j<P-l (4.55) 
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or 


M-l P-1 


RxrXiP^r[n] - s; Mr[n] - Mj[n] + j - s]hl^''[s 


r=0 s=0 

= Rdxi[n;n- Mj[n] + j], 

0<i<M-l,0<j <P-1, (4.56) 

where p = (n) ji_. These are called the discrete multi-channel, multirate Wiener-Hopf 
equations. 

As in the ordinary Wiener hltering problem, the multi-channel, multi¬ 
rate Wiener-Hopf equations can be derived in matrix form by writing (4.54) as 

M-l 


d[n\ = y^(xi [Mj[n]])^hfb 


i=0 


where if 




n = 


Xi[Mi[n] - P -h l] 
Xi[Mi[n] - P] 


X. 




then 




n = 


Xi 

Xi\Mi[n] - l] 


Xi\JAi[n] - P +1] 

with Xj [Mj [n]] as dehned in Section II.B.3, and where is dehned as 

hflO] 

A?’11] 


= 


h^\P-l] 


(4.57) 


(4.58) 


(4.59) 


(4.60) 


where p = (n) ji_. 
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Applying the orthogonality principle yields: 


M-l 


£{xi[Mj[n]]e*[n]} = £ <| | d*[n] - (x^j j ^ ^ 

0<i<M-l, (4.61) 


r=0 


and the Wiener-Hopf equations can be expressed in their matrix form as 

M-l 


^ ^XiXr Ijvt. [„]1 

r=0 


, 0 < i < M - 1, 


(4.62) 


where p = (n) ji_, and 

''.^[jvt.w] =«{“:["K‘[M.N]}. (4.63) 

Applying the second part of the orthogonality principle yields an ex¬ 
pression for the minimum mean-square error: 


. [n 

&min >• -J 


M-l 


E{d[n]e*[n]} = £ ^ d[n] j d*[n] — (x^[Mr[n]] 

0 < i < M - 1, (4.64) 


r=0 


(jp = i?d[n; 0] - ^ [n; n - Mr[n]] 0 < i < M - 1, (4.65) 

r=0 

where p = and its index on the mean-square error reflects the periodically 

time-varying nature of the error. 

d. Matrix Approach to the Multi-channel, Multirate 
Wiener-Hopf Equations 

In the multi-channel, multirate optimal hltering problem, we also de¬ 
velop a matrix-based method to develop the Wiener-Hopf equations. In this problem, 
M observation signals are available to form the desired signal d, where M > 2. The 
estimate d is formed by summing the output of M periodic hlters of order P and is 
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given by 


M-l 

dk[n] = ^ Q < k < K — 1, (4.66) 

i=0 

where K is the period of the hlters required for estimation and is given by = (n)^, 
where N is the system period. For any index n, the observation signal Xi, represented 
by the sequence Xi [n] , can be expressed as 

xf^[n] = D^^Xj[n], 0 </c < iF* - 1, (4.67) 

where 

Xi[n] 

- r 1 Xi[n-l] 

x4nj = 

Xi[n - PKi + 1 ] 

with Xj[n] as dehned in Section II.B.3, and is a decimation matrix with delay 
(3.60), where Ki is the decimation factor associated with the channel, and P is 
the hlter order. Again, note that, in this formulation, all signals and computations 
are at the system rate, and, as a result, we are solving a more general problem since 
we estimate d[n]T^ at every point on not just on thT'd- The desired estimate can 
be recovered by decimating the result by Kd. 

Again, in the usual manner, we form the error in estimation and evoke 
the orthogonality principle, 

d[n] - ^(5ff^[n])^hf^j | = 0 , 

0<i<M-l, (4.68) 
and the Wiener-Hopf equation (4.30) can be expressed in its matrix form as 

M-l 

E 0 < » < M - 1, (4,69) 

j=0 
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and k = {n)j^, or in terms of single-rate statistical parameters (see (4.16) and (4.17)), 


M-l 


hf * = 0 < ^ < M - 1. 

j=0 


(4.70) 


Applying the second part of the orthogonality principle yields an ex¬ 
pression for the minimnm mean-square error given by 

M-l 




^.Jn\ = e{d[n\e*[n\} = £ <j d[n\ d[n\ - ^(5ff^[n])^hf^ 


i=0 


or 


M-l 


a'i = 0] - ht’A!'"*, 0 < » < M - 1, 

i=o 

or in terms of single-rate statistical parameters, 


M-l 


4 = 


Rd4 0] - E 0 < i < M - 1, 

j=0 


(4.71) 


(4.72) 


(4.73) 


where the subscript k on the mean-square error reflects the periodically time-varying 
nature of the error. 


D. CHAPTER SUMMARY 

In this chapter, the problem of signal estimation is introduced from an optimal 
hltering perspective. Specihcally, we consider linear hltering, where the optimality 
criterion is based on minimizing the mean-square error (Wiener hltering). 

To begin the discussion, the Wiener-Hopf equations for single-rate sequences 
are reviewed, establishing the foundation upon which the single- and multi-channel, 
multirate Wiener-Hopf equations are developed. These results are important to the 
methods of signal and image reconstruction discussed in Chapter V. 

Additionally, the concept of index mapping is proposed, which systematically 
describes the relationship between samples of a signal at one rate to those of a signal 
at another rate (indices in one signal domain to those in another signal domain). This 
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concept is particularly well-suited for developing the required regions of support in 
linear filtering. In detailing this relationship, both a causal and non-causal mapping 
transformation are developed, but the notion of mapping is generalized, providing 
a basis for the development of “generalized” multi-channel, multirate Wiener-Hopf 
equations. 

Finally, matrix-based approaches, using the matrix representations of Chap¬ 
ter III, are used to develop the single- and multi-channel, multirate Wiener-Hopf 
equations for implementation in Chapter V. 
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V. SUPER-RESOLUTION SIGNAL AND 
IMAGE RECONSTRUCTION 

In this chapter, we apply the multirate and optimal estimation theory results 
of Chapters III and IV to the problem of signal and image reconstruction. We start 
by developing a 1-D methodology and begin the analysis of this method on a simple 
test signal (a triangle wave) and evaluate the performance on this known data. We 
then apply the procedure to image data when the image is processed along rows. 

We then turn to the full problem of image reconstruction and discuss how 
the 1-D multirate theory and methods extend to 2-D, present an image reconstruc¬ 
tion methodology, and provide an example of the application of this methodology, 
comparing results to other methods. 

A. SIGNAL RECONSTRUCTION 

Consider the problem of estimating a discrete random process d[n], which 
cannot be observed directly, from a set of M related observation signals, represented 
by {a;o[mo], a;i[mi],..., These observation signals are related to the 

random process d[n] through various forms of distortion and interference. Further, 
these signals may be at rates different from that of d[n], and observation signals at 
the same rate may also be shifted in time with respect to one another. 

1. Observation Model 

In order to facilitate discussion, we present the observation model depicted in 
Figure 5.1, which represents the observation signal. In this model, we consider 
linear forms of distortion. Notice that the translation is along the system grid (T^r) 
and is indexed to the particular signal of interest (see e.g., sq and si in Figure 5.2). 
Further, notice that the downsampling factor is also indexed to a particular signal. 
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Distortion 


Translation 


Downsampling 


d{n\ 



TliW 


Additive Noise 


Figure 5.1. Observation model, where observation signals Xi[mi\ are derived from an 
underlying signal d, subject to distortion, additive noise, translation, and downsam¬ 
pling. 

If we consider signal S 3 [n], its associated downsampling factor is L 3 , which may be 
different from that of other observation signals. 

2. Optimal Estimation 

As we have previously observed, if the desired signal d[n\ and its observations 
Xj[mJ are jointly wide-sense stationary, then the linear hlters required for optimal 
mean-square estimation are periodically time-varying [Ref. 7]. In this case, we can 
write the estimate as 

M-l 

dk [n] = (if ^ [n] )^hf ^ (5.1) 

i=0 

where hf^ is a set of time-varying hlter coefficients of length Pi and Xj[n] is a vector 
of samples from the observation sequence. The periodic time variation is denoted 
by the index k where 0 <A;<L — l,Lis the system periodicity and k = {n) ^ (see 
Figure 5.3). In this analysis, we consider that the observations signals are maximally 
decimated (see Section III.C.4) versions of d[n], where Li = L, and, in this case, the 
observation vectors can be written as 

(5.2) 
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nTj 



Figure 5.2. Observation sequences sq and si shifted by a delay (f = 0, i = 1, respec¬ 
tively). 


and 

Sj[n] = [si[n],Si[n - 1],.. .,Si[n - PiL + 1]]^ . (5.3) 

The decimation matrix with time delay extracts the appropriate samples 

from Sj[n] to form each observation vector. 

Minimizing the mean-sqnare error [Ref. 9] leads to a set of Wiener-Hopf 
equations of the form 


■p (^) 

-■^oo 

p (^) 

Kqi 

1 

7 

o 


1 

* 


~ik)* 

'■do 

p {k)*T 

p (^) 

Kii 

p (^) 


hW* 

= 

~ik)* 

'-dl 

p(fc)*r 

p(fc)*r 

p (k) 

■ ^L-1L-1_ 


1 (k)* 

[hi-ij 


~ik)* 

_ dL-l_ 


0<fc<L-l, (5.4) 

where the time average mean-sqnare error [Ref. 9] is given by 

( 5 - 5 ) 

fc=0 j=0 
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v(‘) 


[«] 


xf\n] 



k =(n) 


Figure 5.3. Reconstruction of the original signal from an ensemble of subsampled 
signals based on optimal linear filtering 


Notice that this is an arithmetic mean of the set of al, 0 < k < L — 1. The correlation 
terms are dehned as 

-IX, _ (5.6) 

rW = (5.7) 

and 

Rd[0] = £,{d[n]d*[n]} (5.8) 

where 

Vdi = £{d[n]s. [n]} (5.9) 

and 

Rij = £{s4n]s*^[n]}. (5.10) 

Solving the multirate Weiner-Hopf equations (5.4) yields a set of hlter coefficients, 
which can be used in the estimation of d[n] as depicted in Figure 5.3. The application 
of these hlters to the observation data is illustrated in Figure 5.4. 
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Figure 5.4. Reconstruction of the original signal from an ensemble of subsampled 
signals based on FIR Weiner filtering with decimation factor L = 3 and hlter order P 
= 4. The hgure illustrates the support of the time-varying hlters at a particular 
time, n = 15 and fc = 0 (shaded circle). 


The system equations can alternatively be derived using the appropriate map¬ 
ping function. In this case, we hnd that the causal mapping function can be written 


M=r : n 


n — i 


n 

+ 

k — i 

L 


[l\ 

L 


( 6 . 11 ) 


or 


mi = Msr(n), (5.12) 

where L is the decimation factor, and k = {n )Notice that the hrst term of 
the mapping function corresponds to the multirate hltering scenarios of Chapter IV; 
however, we require an additional corrective term to account for the translation 
between signals. 

Notice that if 


n — i 


n 

+ 

k — i 

L 


[l\ 

L 


using the dehnition (2.41) of the common residue (n)^^, we can write 


n = 


n 
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Then, by direct substitution 


n — i 


n i 


[f J L in), ^ 

L 


_L ~ L_ 


1 

b^ 



Since |_^J is always an integer, it can be moved out of the floor operation (2.40), to 
obtain 


n — i 


n 

+ 

Wl 

i 

L 


-L. 

L 

~ L_ 


By dehnition, k = {n)j^] therefore, we can write 


n — i 


n 

+ 

k — i 

L 


.L. 

L 


Example 13. As an example, for L = 3, we have the following mappings from a 
partieular estimate index to the corresponding observation indices {recall that for the 
maximally decimated case with L = 3, i = {0, 1,2}) as shown in Table 5.1. 


n 

mi = Msr(n) 

12 


13 

1 |co 

+ 

14 

{[tJ + [¥J} = {3.4,4} 


Table 5.1. Causal mapping from an estimate signal index to the associated observation 
signal indices, for the maximally-decimated case, L=3. 


These mapping transformations can be directly applied to the multi-channel, multi¬ 
rate Wiener-Hopf equations (4.62) and (4.65). 

3. Reconstruction Methodology 

The process used in the reconstruction of a signal from a set of subsampled 
observations is described as follows. First, a training signal, at the desired rate, is 
obtained that is representative of the class of signals that will be processed. From 
this signal, a maximally-decimated set of observations are derived and the single¬ 
rate statistical parameters are extracted. Then, using the Wiener-Hopf equations 
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developed in Section V.A.2, the filter coefficients are computed. With the class- 
specific filter coefficients, we are able to reconstruct signals “of the same class” by 
employing the estimate of (5.1). In other words, we use training data to develop 
the required filter coefficients from a class-representative signal at the desired rate, 
and apply these filters to any set of “representative” observations at a lower rate to 
reconstruct the desired signal. 

4. Application results 

a. Reconstruction of a Known Waveform 

To evaluate the performance of the proposed method, two examples are 
presented. In the first example, a triangular waveform is considered for reconstruction. 
Our method was compared to the method described in [Ref. 31], which can produce 
an exact reconstruction of the triangular waveform if the highest frequency terms are 
left out. Both methods produce accurate results when there is no noise added to the 
observation sequences. When a small amount of noise is added to the observation 
sequences, the exact reconstruction method fails to reliably reproduce the signal while 
the method described here continues to produce a reasonably good approximation to 
the signal even under severely noisy conditions. The original triangular waveform is 
shown in Figure 5.5 (a). Also shown there are the reconstructed waveform (a) and 
the mean-square error (b). 

The observation sequences are shown in Figure 5.6. These were con¬ 
structed by shifting the original sequence, downsampling by a factor of L = 3, and 
adding white Gaussian noise. In this particular example, a signal-to-noise ratio (SNR) 
of -4.8dB was used. For our purposes, SNR is computed from the ratio of signal power 
to noise variance. Note that the underlying form of the original sequence is unde¬ 
tectable. 
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SNR= -4.8dB, Filter Length P= 8, Decimation Factor Z, = 3 




Figure 5.5. Simulation results using optimal linear filtering method for reconstruction; 
SNR = —4.8dB, P = 8, and L = 3. 





Figure 5.6. Observation sequences of an underlying triangle waveform after being 
subjected to additive white gaussian noise and subsampled by a factor of L = 3. 
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As a precursor to consideration of the two-dimensional image recon¬ 
struction problem, we applied the one-dimensional methods described here to image 
data. 

b. Extension to Two-Dimensional Reconstruction 

We consider each row of the observed LR images as an observation sig¬ 
nal vector belonging to the set {xo[m], xi[m],... ,xm-i[?7i]} (see Figure 5.7). Recon¬ 
struction is then accomplished line-by-line until every row of each image is processed. 
In this case, the original image is depicted in the left panel of Figure 5.8, and one of 



Figure 5.7. Line-by-line processing of observation images. 

its three subsampled observation images with additive white Gaussian noise is given 
in the right panel. The image depicted in the left panel of Figure 5.9 represents 
the result of applying standard nearest-neighbor interpolation to one of the three 
noisy subsampled images, and the reconstructed image is shown on the right panel. 
Although the image is processed in only one direction, there is signihcant improve¬ 
ment over the interpolated image. In particular, note that edges of structures can be 
observed in many cases where the interpolated image does not provide such detail. 
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Figure 5.8. Original image (left); Subsampled observation image L = 3 with additive 
white Gaussian noise, OdB (right). 



Figure 5.9. Result of applying standard nearest-neighbor interpolation to one of the 
three noisy subsampled images (left); Reconstructed image (right). 
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B. IMAGE RECONSTRUCTION 


In this section, we consider the signal reconstruction problem in two dimen¬ 
sions, i.e., super-resolution (SR) image reconstruction. This analysis is an extension 
of the theory and application presented for the 1-D case. Specifically, we consider 
the problem of estimating a two-dimensional discrete random process, which cannot 
be observed directly, from a set of related observations, at a different sampling rate. 
Like the 1-D problem posed earlier, the observation signals are related to the random 
process through various forms of distortion and interference and these signals may 
also be shifted in time with respect to one another. 

1. Observation Model 

In the context of SR image reconstruction, the underlying two-dimensional 
signal (image) is at a high-rate (HR), and the observations are at some low-rate 
(LR). The relationship between the HR signal and a particular observation image 
((b j)*'^ ~ channel) is illustrated in the block diagram of Figure 5.10. This observation 
model shows that each LR observation is acquired from the HR image subject to 
distortion (typically blur), subpixel translation, downsampling, and channel noise. 
We can represent the observation model as 

F., = DgFG.pg’' + (5.13) 

The matrix F G represents the desired HR two-dimensional signal (image) 

with {MNL 1 L 2 ) pixels. Its related LR observation is represented by the matrix 
Fjj G The matrix Gij is a linear operator that accounts for channel distortion 

and matrix represents the channel noise. The parameters Li and L 2 represent the 
horizontal and vertical downsampling factors, respectively, and the parameters i and 
j represent the horizontal and vertical subpixel translation, respectively. The matrix 
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F 



Fij 


t 

mj 


Figure 5.10. Observation model relating the HR image with an associated LR ob¬ 
servation. Each LR observation is acquired from the HR image subject to distortion 
(typically blur), subpixel translation, downsampling, and channel noise. 

is a decimation matrix with delay, dehned in (3.60), and is used to extract the 
appropriate pixels from the HR image to form a particular observation image. The 
matrix on the left incorporates the horizontal downsampling and translation (Li,z) 
while the matrix on the right incorporates the vertical downsampling and translation 
(L 2 , j)- A maximally decimated set (see Section HI.C.4) of images is obtained from 

{Fy }, f = 0,1,..., Li - 1, j = 0,1,..., La - 1. 


Example 14. This example illustrates the operation of decimation matrices with 
delay on a HR image. Consider Li = L 2 = 3, i = 1 and j = 2; the HR image F with 
Gi ,2 = I and r]i 2 = 0 : 
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The product Fi ^2 can he written 
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The observation imageY 1^2 is one of nine (L 1 XL 2 ) possible configurations. ■ 

2. Optimal Estimation 

Consider the set of LR observations {Fjj} and the related HR image F. We 
desire to form an estimate for the HR image by some weighted sum of the LR obser¬ 
vations. The estimate can be written as 

F[ni,n2] = ^ (5.14) 

i=0 j=0 

where the expression on the right represents the Frobenius inner product of the ma¬ 
trices , dehned in Dehnition 12. Here, the matrix tij[ni,n 2 ] G is the set of pixels 

that form the region of support within the appropriate LR observation, required to 
form this estimate, and matrix is the corresponding set of Liter coefficients. In 

this discussion, this region of support is called the image mask and its associated set 
of filter coefficients is called the filter mask. The filter masks are chosen to minimize 
the mean-square error 

£{||F-Ff}, (5.15) 

where the norm is the Frobenius norm; £{■} is the expectation operator; ki = 
for 0 < ki < Li] and {i,j) represent subpixel translation, with 0 < i < Li — 1, 

0 < j < L 2 — 1, and {i,j) G Z. In the maximally-decimated case, i and j span the 
entire set {0,1,..., Li — 1} and {0,1,..., L 2 — 1}, respectively. Both the set of LR 
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image masks{fjj[ni, 77 , 2 ]} and filter masks are further described in Sections 

V.B.2(a) and V.B.2(c), respectively. 

a. Index Mapping 

Developing the LR image masks {fp} involves mapping indices in the 
HR sample index domain d'i? to those in the LR sample index domain d'p. For a 
given pixel intensity F[ni,n 2 ], the HR indices [ni,n 2 ] map to a set of LR sample 
indices {[mi,m 2 ]p}, which correspond to pixel intensities {^^[ 7711 ,m 2 ]}. The indices 
corresponding to each observation matrix Fp are determined such that 

|[77i,772] - [mi,m 2 ]p| (5.16) 

is minimized for each 

This mapping can be shown to be 


[mi,m2]p = [M*(77 i),Mj(772)], 


(5.17) 


where M[77] is dehned 


n — i+ |_-|J 


n 

+ 

k — i+ |_-|J 

L 


-L. 


L 


mi = Msr(77), 


(5.18) 

(5.19) 


where nTd E m^Tj G Tp k = (n) and L is the decimation factor. Notice that the 
hrst term of the mapping function corresponds to the multirate hltering scenarios 
of Chapter IV; however, we require an additional corrective term |_^J to account 
for the translation between signals. 

b. LR Image Mask 

The LR indices [mi,m 2 ]p for each observation represent the centroid 
of each of the LR image masks. Given a desired mask size of P x Q, each LR image 
mask is comprised of the P x Q pixels closest to Fp[mi,m 2 ]. 

c. Filter Mask 

If the desired HR image F and its observations {Fp} are jointly homo¬ 
geneous, then the linear Liters required for optimal estimation are periodically 
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spatially-varying, an extension of [Ref. 7]. This periodicity can be described in terms 
of the “phase” (/ci,/c 2 ), where ki = njmodLj. If we dehne the set of least posi¬ 
tive residues as = 0,1 ,..., — 1, then ki G and k 2 G A^j, and all possible 

combinations of phase can be represented as A^^ x 

Figures (5.11) and (5.12) depict the phase variation for Li = L 2 = 2. 
In this case, A 2 = {0,1} and A 2 x A 2 = {(0, 0), (0,1), (1, 0), (1,1)}. The spatial 
periodicity of the phase can be observed by noting the regular recurrence of phase 
terms. 
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Figure 5.11. Index representation to modulo representation with Li = L 2 = 2 (note 
the spatial phase periodicity). 


3. Reconstruction Methodology 

a. Least Squares Formulation 

In order to determine the hlter masks required to estimate 

the HR image, a least squares (LS) approach is employed. We identify the set of all 
HR pixels that correspond to a given phase and denote this set of HR pixels as the 
matrix This concept is depicted in Figure 5.12 where each shape corresponds 

to a unique phase (circle (0,0), square (0,1), triangle (1,0), and star (1,1)). From 
(5.14), we can see that 


F[/i,/2] 

F[mi,m2] 

F[ni,n2] 


h' 

ij 


(5.20) 
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Figure 5.12. Relationship between HR pixels and spatially-varying filter masks in 
formulating the LS problem with Li = L 2 = 2. 


where the left hand expression is From this, we can write 

YikiM) (5.21) 

where $ is the data matrix. This system of equations is solved in a least squares 
sense for the required set of hlter masks at each phase 

b. Processing Method 

The process used in the SR image reconstruction of a set of LR ob¬ 
servations is described as follows. First, a HR training image is obtained that is 
representative of the class of images that will be processed. From this image, a 
maximally-decimated set of LR observations are derived and then through the least 
squares methodology of Section V.B.3(a), hlter coefficients are computed. With the 
class-specihc hlter coefficients, we are able to reconstruct images “of the same class” 
by employing the estimate of (5.14). In other words, we use training data to develop 
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filter masks from a class-representative HR image, and then apply these filters to any 
set of “representative” LR observations to reconstrnct a HR image. 

4. Application Results 

In order to evalnate the performance of this method of image reconstrnction, 
we process the “skyline” image depicted in Fignre 5.13, snbject to varying degrees of 
AWGN. The image nsed for the training process is the 204 x 204 pixel image segment 
depicted in this figure. From this image, a set of filter coefficients is derived that is 
used in SR image reconstruction. 



Figure 5.13. Image segment used to train filter. 

The target or object of the reconstruction is depicted in Figure 5.14. From 
this 204 X 204 pixel image segment, LR observations are derived and are hltered using 
the class-specihc hlter masks. The same level of AWGN is used for the training and 
the target images. 

Figure 5.15 depicts three members of the set of LR observations with various 
subpixel translations. The first image represents subpixel translation by one pixel 
in the horizontal direction and no translation in the vertical direction. The second 
represents translation by one pixel in both the vertical and horizontal directions. The 
third image represents translation by two pixels in both directions. 
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Figure 5.14. Image segment to be estimated. 


In the remaining figures, the left panel depicts the SR image reconstruction 
using the proposed algorithm and the right panel depicts nearest-neighbor interpola¬ 
tion of one of the LR observations. In every case, the proposed method is superior to 
the interpolated result. During these experiments, other interpolation methods were 
considered, including bilinear and bicubic methods; again, the proposed method was 
found to be superior. 

Figure 5.16 compares reconstructed and interpolated images for the case of no 
additive noise, downsampling by 3 in both the vertical and horizontal directions, and 
with hlter mask size of 3 x 3. In this case, the reconstruction yields a result that is 
visibly indiscernible from the target image. 

Figure 5.17 compares reconstructed and interpolated images for the case of an 
SNR = 5 dB, downsampling by 3 in both the vertical and horizontal directions, and 
with filter mask size of 3 x 3. In this case, we see the effects of additive noise on the 
reconstruction. Despite some blurring of edges, details are still discernible. 

Figure 5.18 compares reconstructed and interpolated images for the case of an 
SNR = —1.5 dB, downsampling by 3 in both the vertical and horizontal directions, 
and with filter mask size of 3 x 3. In this case, the effects of additive noise on 
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Figure 5.15. Downsampled observation images with subpixel translations (1, 0), (1,1), 
and (2, 2), respectively; Li = L 2 = 3, P = Q = 3, and no AWGN. 

the reconstruction are quite deleterious. Further blurring of edges is evident and 
details have become hard to see. However, major features in the image are still 
discernible. Thus, there is a significant advantage in using the proposed method 
over interpolation, where not even major features are discernible. Note that the 
poorer performance of the interpolation methods is not unexpected. In this case, 
only a single LR observation image is used for reconstruction while in the case of 
the proposed method, multiple, independent LR observations are used. Despite this 
obvious disadvantage, the interpolation methods are used throughout the literature 
as a benchmark [Ref. 1, 2, 28] and we follow here. 

C. CHAPTER SUMMARY 

In this chapter, the signal and image reconstruction problem is considered 
from a multirate point of view (Chapter III), using the multirate optimal estimation 
theory presented in Chapter IV. 

First, the problem of reconstruction is posed in one dimension, in terms of 
a set of observation sequences that are related to a desired random sequence. An 
observation model is presented that models this relationship, accounting for linear 
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Figure 5.16. Comparison between a reconstructed image and interpolated image; 
Li = La = 3, P = Q = 3, no AWGN. 

distortion, additive noise, downsampling, and subsample translation. Signal recon¬ 
struction is achieved by forming an estimate comprised of a linear combination of 
samples of the observation sequences, and then developing and solving a set of re¬ 
lated linear equations for the required hlter coefficients. These coefficients are used 
to hlter LR data. Finally, the results of this proposed methodology are presented and 
compared to other reconstruction methods. 

Next, the super-resolution reconstruction problem is considered and posed in 
terms of a set of 2-D observation sequences that are related to a desired 2-D random 
sequence. Again, an observation model is presented that models this relationship, 
accounting for linear distortion, additive noise, downsampling, and subpixel transla¬ 
tion. Image reconstruction is achieved by forming a linear estimate comprised of a 
linear combination of pixels in each LR observation image mask and developing and 
solving a set of related linear equations for the required hlter coefficients. Formation 
of these image masks is based on extension of the non-causal 1-D index mapping 
results to two dimensions. The resultant hlter masks are used to hlter LR image 
data. Finally, the results of the proposed methodology are presented and compared 
to other reconstruction methods. 
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Figure 5.17. Comparison between a reconstructed image and interpolated image; 
Li = L 2 = 3, P = Q = 3, and SNR = 5 dB. 



Figure 5.18. Comparison between a reconstructed image and interpolated image; 
Li = L 2 = 3, P = Q = 3, and SNR = -1.5 dB. 
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VI. CONCLUSION AND FUTURE WORK 


A. SUMMARY 

In this dissertation, new signal processing methods for multirate signals in one 
and two dimensions are developed and applied to the problem of “super-resolution” 
signal and image reconstruction. In super-resolution processing a signal or image with 
high-resolution is constructed by using data from several images of a lower resolution. 

A signihcant contribution of this work is the development of the theory of mul¬ 
tirate systems, which provides the foundation upon which the proposed reconstruction 
methods are built. In developing this theory, a conceptual framework for the anal¬ 
ysis of multirate systems is cited, which enables the extension of optimal estimation 
and linear hltering theory to multirate systems. Further, it leads to the concept of a 
system grid, which allows representation of the various signals in a multirate system 
on a common domain. 

System theory for multirate systems, also developed here introduces multirate 
adaptations of the concepts of shift-invariance, periodic shift-invariance and causality. 
In addition, the input-output relation of a multirate system is discussed in terms of the 
system response and its associated Green’s function and is adapted to the multirate 
problem in both system-level and signal-level representations. Finally, a method 
is adopted to provide matrix representations for the operations of downsampling, 
upsampling and linear hltering. 

Another signihcant contibution of this work is the treatment of the problem 
of signal reconstruction from a multirate, optimal hltering perspective. Recognition 
that signal or image reconstruction can be viewed as a problem in signal estimation, 
where a related low-rate (low-resolution) signal or set of signals is used to estimate 
an underlying high-rate (high-resolution) signal, and that the observation signal or 
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signals and the desired signal form a multirate system, motivates this work. In partic¬ 
ular, it suggests that extension of well-known results for single-rate systems could be 
extended to multirate signals and systems. In development of the resulting multirate 
forms of the Wiener-Hopf equations, the concept of index mapping is proposed. Index 
mapping systematically describes the relationship between samples of a signal at one 
rate to those of a signal at another rate and is particularly important for developing 
the required regions of support in linear hltering. In developing this relationship, both 
causal and non-causal mapping transformations are considered. The overall concept 
is more general, however, and provides a basis for the development of “generalized” 
multi-channel, multirate Wiener-Hopf equations. 

In applying the multirate and optimal estimation theory results to the prob¬ 
lem of signal and image reconstruction, a one-dimensional method is hrst explored. 
The analysis of this method starts with a simple test signal (a triangle wave) and 
its performance is evaluated on this known data. The procedure is then applied to 
image data when the image is processed along rows. Finally, the full two-dimensional 
problem of image reconstruction is considered and an image reconstruction method¬ 
ology is developed. An example of the application of this method is provided with a 
comparison to results of other methods. 

B. FUTURE WORK 

As fundamental limits are reached in future image acquisition systems, inter¬ 
est in signal processing solutions will continue to intensify, and in particular, the area 
of super-resolution image reconstruction will likely be of great interest. During this 
research, advances were made in the areas of multirate signal processing and optimal 
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estimation theory, which led to development of the proposed reconstruction methods. 
With these hndings, other research opportunities became apparent, and are mentioned 
here. 

In the development of the multirate signal processing theory, several num¬ 
ber theoretic results were obtained to describe the relationship between constituent 
signals in a multirate system. These results were shown to extend directly to the 
second moment analysis of such systems, but the second moment relationships were 
never entirely understood. Future work in characterizing second moments, and even 
higher moments, would provide further insight into the behavior of multirate systems. 
Furthermore, frequency-domain characterizations of such systems would also provide 
valuable insight, and would be a useful extension of [Ref. 12] on bifrequency and 
bispectrum maps. 

In the analysis of performance, the 2-D method developed here was compared 
to standard interpolation methods, using a particular class of images. In future 
work, it would be benehcial to compare the proposed method to other “state of 
the art” algorithms [Ref. 1, 2]. Further, a more robust analysis of the proposed 
method’s performance would be achieved if several classes of images were available 
for comparison. Finally, it would be useful to examine the effects of linear distortion, 
image rotation and non-Gaussian noise on performance. 

Finally, in conducting this research, limitations were imposed on the sampling 
rate of constituent signals of a multirate system. For our purposes, sampling rates 
were constrained to be integer-valued. The theory developed in this work and in 
[Ref. 59] could be generalized to include cases where real-valued sampling rates are 
considered. 
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